Arts, Faculty of
Geography, Department of
DSpace
UBCV
Ng, Darian Chi Yung
Master of Science - MSc
University of British Columbia
Freshwater wetlands are effective nature-based solutions to climate change for their unique ability to sequester atmospheric carbon for extended periods of time. However, they simultaneously emit large amounts of methane gas (CH₄). CH₄ is a highly potent greenhouse gas, and its spatial heterogeneity in wetlands is not yet fully understood. This study is the first of its kind and explores a novel approach to characterizing and predicting the spatial distribution of CH₄ fluxes (FCH₄) within a single flux tower footprint using a combination of remote-sensing via the Landsat 8 satellite and eddy covariance (EC) Footprint-weighted Flux Maps. There were two research objectives: 1) to determine whether there exist significant relationships correlating remote-sensing with within-footprint FCH₄, and 2) to investigate the potential of extrapolating and upscaling these relationships to predict FCH₄ beyond the flux tower. The Normalized Difference Vegetation Index (NDVI), Water Index (NDWI), Moisture Index (NDMI), and Land Surface Temperature (LST) remote-sensing products were collected from Landsat 8, and EC data from four wetland sites, MBPPW1, MBPPW2, US-Myb, and US-WPT were used for the Flux Maps. Hierarchical clustering was used on the Flux Maps to identify distinct FCH₄ features, and then applied to Landsat 8 to process the data for use in two stages of regression analysis. The first stage used simple linear regression to assess the relationships between FCH₄ and the four remote-sensing indices. The second stage used multiple linear regression (MLR) to model the spatial distribution of FCH₄. The simple linear regression analysis demonstrated strong relationships between remotely-sensed surface conditions and FCH₄, and found stronger FCH₄ from open water at MBPPW2 and US-WPT, and stronger FCH₄ from the emergent vegetation at US-Myb. MBPPW1 demonstrated a limitation of this approach requiring sites to have strong FCH₄ signals. The MLR models supported the feasibility of upscaling using remote-sensing by capturing much of the FCH₄ variance, but had reduced accuracy from strong error-biases and low regression slopes. This study demonstrates the applicability of its novel methodology to upscaling FCH₄ using remotely-sensed data and paves the way for future studies aiming to improve FCH₄ estimates from wetlands at a regional level.
CHARACTERIZATION AND PREDICTION OF METHANE SPATIAL HETEROGENEITY WITHIN AND BEYOND A FLUX TOWER FOOTPRINT by Darian Chi Yung Ng B.Sc., The University of British Columbia A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Geography) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2023 © Darian Chi Yung Ng, 2023 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled: Characterization and prediction of methane spatial heterogeneity within and beyond the flux tower footprint. Submitted by Darian Chi Yung Ng in partial fulfillment of the requirements for the degree of Master of Science in Geography Examining Committee: Dr. Sara Knox, Assistant Professor, Geography, UBC Co-supervisor Dr. T. Andrew Black, Professor, Land and Food Systems, UBC Co-supervisor Dr. Camilo Rey-Sanchez, Assistant Professor, Marine Earth and Atmospheric Sciences, North Carolina State University Supervisory Committee Member Dr. Brett Eaton, Professor, Geography, UBC Additional Examiner iii Abstract Freshwater wetlands are effective nature-based solutions to climate change for their unique ability to sequester atmospheric carbon for extended periods of time. However, they simultaneously emit large amounts of methane gas (CH4). CH4 is a highly potent greenhouse gas, and its spatial heterogeneity in wetlands is not yet fully understood. This study is the first of its kind and explores a novel approach to characterizing and predicting the spatial distribution of CH4 fluxes (FCH4) within a single flux tower footprint using a combination of remote-sensing via the Landsat 8 satellite and eddy covariance (EC) Footprint-weighted Flux Maps. There were two research objectives: 1) to determine whether there exist significant relationships correlating remote-sensing with within-footprint FCH4, and 2) to investigate the potential of extrapolating and upscaling these relationships to predict FCH4 beyond the flux tower. The Normalized Difference Vegetation Index (NDVI), Water Index (NDWI), Moisture Index (NDMI), and Land Surface Temperature (LST) remote-sensing products were collected from Landsat 8, and EC data from four wetland sites, MBPPW1, MBPPW2, US-Myb, and US-WPT were used for the Flux Maps. Hierarchical clustering was used on the Flux Maps to identify distinct FCH4 features, and then applied to Landsat 8 to process the data for use in two stages of regression analysis. The first stage used simple linear regression to assess the relationships between FCH4 and the four remote-sensing indices. The second stage used multiple linear regression (MLR) to model the spatial distribution of FCH4. The simple linear regression analysis demonstrated strong relationships between remotely-sensed surface conditions and FCH4, and found stronger FCH4 from open water at MBPPW2 and US-WPT, and stronger FCH4 from the emergent vegetation at US-Myb. MBPPW1 demonstrated a limitation of this approach requiring sites to have strong FCH4 signals. The MLR models supported the feasibility of upscaling using remote-sensing by iv capturing much of the FCH4 variance, but had reduced accuracy from strong error-biases and low regression slopes. This study demonstrates the applicability of its novel methodology to upscaling FCH4 using remotely-sensed data and paves the way for future studies aiming to improve FCH4 estimates from wetlands at a regional level. v Lay Summary Freshwater wetlands are capable of removing large amounts of carbon from the atmosphere and storing them for long periods of time in their flooded soils. At the same time, however, they are also strong producers of methane gas (CH4). CH4 is a more potent greenhouse gas than carbon dioxide, so it is important to understand CH4 fluxes (FCH4) from wetlands, but the spatial patterns of their emissions are not yet fully understood. This thesis presents a novel approach to addressing this knowledge gap by combining satellite data with surface-flux maps to spatially characterize FCH4. The results were able to strongly correlate satellite data with FCH4 at four studied wetland sites. From here, multilinear regression models were able to predict FCH4 with moderate accuracy. These findings form a strong framework for the development of a valuable tool capable of leveraging open-access satellite data to improve estimates of FCH4 in freshwater wetlands. vi Preface I was responsible for processing the high frequency flux data for MBPPW1 and MBPPW2. The code for this was provided by my supervisor, Dr. Sara Knox, and the REddyProc (v2.2.0) R package by Thomas Wutzler. I was greatly assisted by Zoran Nesic, who helped monitor and troubleshoot the data coming from these two sites. Site maintenance and field work at MBPPW1&2 were performed by the team at Ducks Unlimited Canada, led by Dr. Pascal Badiou, Paige Kowal, and Bryan Page. The code for the Footprint-weighted Flux Maps was developed and provided by Dr. Camilo Rey-Sanchez. vii Table of Contents Abstract ......................................................................................................................................... iii Lay Summary .................................................................................................................................v Preface ........................................................................................................................................... vi Table of Contents ........................................................................................................................ vii List of Tables ................................................................................................................................ ix List of Figures .................................................................................................................................x List of Abbreviations .................................................................................................................. xii List of Symbols ........................................................................................................................... xiv Acknowledgements ......................................................................................................................xv Dedication ................................................................................................................................... xvi Chapter 1: Introduction ................................................................................................................1 1.1 Wetland Carbon Cycle .................................................................................................................... 2 1.2 Wetland Methane Emissions .......................................................................................................... 4 1.3 Research objectives .......................................................................................................................... 6 Chapter 2: Methods .....................................................................................................................10 2.1 Study sites ............................................................................................................................10 2.2 Eddy covariance data ...........................................................................................................14 2.2.1 Eddy covariance measurements at MBPPW1 and MBPPW2 ......................................15 2.2.1.1 Data filtering ..............................................................................................................16 2.2.1.2 Data gap-filling ..........................................................................................................17 2.2.1.2a Gap-filling meteorological data ...............................................................................17 2.2.1.2b Gap-filling CO2 and heat flux data ..........................................................................18 2.2.1.2c Gap filling FCH4......................................................................................................18 2.3 Footprint-weighted Flux Maps ............................................................................................20 2.4 Landsat 8 Remote Sensing ...................................................................................................22 2.5 Map data processing ............................................................................................................25 2.5.1 Matching spatial layout and coverage...........................................................................26 2.5.2 Matching pixel resolution .............................................................................................26 2.5.3 Hierarchical clustering ..................................................................................................27 2.6 Linear regression analysis ....................................................................................................30 2.7 Multiple linear regression modelling ...................................................................................31 Chapter 3: Results........................................................................................................................33 3.1.1 Simple linear regression – MBPPW1 ...........................................................................33 3.1.2 Simple linear regression – MBPPW2 ...........................................................................36 3.1.3 Simple linear regression – US-Myb ..............................................................................39 3.1.4 Simple linear regression – US-WPT .............................................................................43 3.1.5 Summary of simple linear regression results ................................................................45 viii 3.2 FCH4 Prediction with Multiple Linear Regression .............................................................47 Chapter 4: Discussion ..................................................................................................................52 4.1.1 Methane fluxes at MBPPW1 ........................................................................................52 4.1.2 Methane fluxes at MBPPW2 ........................................................................................55 4.1.3 Methane fluxes at US-Myb ...........................................................................................59 4.1.4 Methane fluxes at US-WPT .........................................................................................62 4.1.5 Linear regression analysis – key findings .....................................................................63 4.1.5a Linear regression analysis - statistical limitations and future work ............................63 4.1.5b Linear regression analysis – practical limitations and future work ............................64 4.2 MLR modelling ....................................................................................................................66 4.2.1 Applicability of MLR predictions .................................................................................67 4.2.2 Limitations of MLR - temporal dependency ................................................................70 4.2.3 Limitations of MLR - multicollinearity and spatial dependency ..................................71 4.2.4 Additional modelling options .......................................................................................74 Chapter 5: Conclusion .................................................................................................................75 5.1 Research objectives and next steps ..............................................................................................76 References .....................................................................................................................................78 Appendix A – Meteorological variable gap filling ....................................................................87 Appendix B – Additional figures ................................................................................................92 ix List of Tables Table 1 Summary table of CH4 flux drivers and the associated remote sensing products… 23 Table 2 MBPPW1 r-value data table………………………………………………………. 36 Table 3 MBPPW2 r-value data table………………………………………………………. 39 Table 4 US-Myb r-value data table………………………………………………………… 23 Table 5 US-WPT r-value data table………………………………………………………... 45 Table 6 MLR model coefficients…………………………………………………………... 49 Table 7 MLR model prediction results…………………………………………………….. 51 Appendix Table 1 Root mean squared error (RMSE) analysis summary………………….. 89 Appendix Table 2 US-Myb r-value data table at  = 0.15 significance level……………... 94 x List of Figures Figure 1 MBPPW1 Photograph and satellite image……………………………………………. 11 Figure 2 MBPPW2 Photograph and satellite image……………………………………………. 12 Figure 3 US-Myb Photograph and satellite image……………………………………………... 12 Figure 4 US-WPT Photograph and satellite image…………………………………………….. 13 Figure 5 Methane flux time series…………………………………………………………….. 14 Figure 6 (MBPPW1) Random forest driver importance graph………………………………….. 19 Figure 7 (MBPPW2) Random forest driver importance graph………………………………….. 20 Figure 8 Footprint-matched Landsat 8 map…………………………………………………… 26 Figure 9 Example Flux Map resolution coarsening…………………………………………….. 27 Figure 10 Hierarchical clustering dendrogram analysis………………………………………… 29 Figure 11 MBPPW1 simple linear regression…………………………………………………. 35 Figure 12 Boxplots of statistically significant Pearson correlation coefficients for MBPPW1…….36 Figure 13 MBPPW2 simple linear regression…………………………………………………. 38 Figure 14 Boxplots of significant Pearson r-values for MBPPW2……………………………… 39 Figure 15 US-Myb simple linear regression…………………………………………………… 41 Figure 16 Boxplots of statistically significant Pearson correlation coefficients for US-Myb……...42 Figure 17 US-WPT simple linear regression…………………………………………………... 44 Figure 18 Boxplots of statistically significant Pearson correlation coefficients for US-WPT…….. 45 Figure 19 Median r-value heat map…………………………………………………………… 46 Figure 20 Predicted-y and residual plots………………………………………………………. 48 Figure 21 Predicted FCH4 plots………………………………………………………………. 50 Figure 22 Box plot of FCH4 at MBPPW1…………………………………………………….. 54 xi Figure 23 Drone photography of MBPPW2…………………………………………………… 57 Figure 24 Tovi Flux Footprint Allocation (FFA) results for FCH4 at MBPPW2………………… 58 Figure 25 Tovi FFA results for CO2 at MBPPW2……………………………………………... 59 Figure 26 Wind direction plot………………………………………………………………… 61 Figure 27 Footprint-weighted Flux Map output………………………………………………... 62 Figure 28 NDMI vs, NDWI scatter plot……………………………………………………….. 62 Figure 29 MLR prediction on un-clustered data……………………………………………….. 68 Figure 30 MBPPW2 Landsat 8 NDVI images…………………………………………………. 69 Figure 31 MLR prediction results on NDVI-clustered data…………………………………….. 70 Figure 32 MBPPW2 Scatterplot matrix of MLR variables……………………………………... 73 Figure 33 Predicted-y plot with weighted variance model ……………………………………... 74 Appendix Figure 1 Pearson correlation tests for MBPPW1 and Hamiota Station………………... 87 Appendix Figure 2 Pearson correlation tests for MBPPW1 and Oakburn Station………………... 88 Appendix Figure 3 Pearson correlation tests for MBPPW2 and Newdale Station……………….. 90 Appendix Figure 4 Simple linear regression scatterplots for MBPPW1, 2022…………………... 92 Appendix Figure 5 Simple linear regression scatterplots for MBPPW1, 2022…………………... 93 Appendix Figure 6 Pearson r-value boxplots for US-Myb at the α = 0.15 significance level……… 94 Appendix Figure 7 US-Myb (2014) linear regression results………………………………………. 95 xii List of Abbreviations Abbreviations Definition CH4 Methane CO2 Carbon dioxide EC Eddy covariance FCH4 Methane flux GPP Gross Primary Productivity H Sensible heat flux IRGA Infrared gas analyzer LE Latent heat flux LST Land surface Temperature MDS Marginal distribution sampling MLR Multiple linear regression NDVI Normalized difference vegetation index NDMI Normalized difference moisture index NDWI Normalized difference water index NEE Net ecosystem exchange R2 Coefficient of determination RH Relative humidity ROI Region of interest SWIN Incoming shortwave radiation TA Air temperature TOA Top of atmosphere xiii Abbreviations Definition TS Soil temperature WD Wind direction WS Wind speed WTD Water table depth xiv List of Symbols Symbols Definition Ai; Bi; Ci Surface radiative transfer coefficients u horizontal wind speed v lateral wind speed w vertical wind speed x Reflectance in wavelength band x s Surface emissivity X Matrix of predictor variables β Predictor coefficient vector Z Matrix of interacting predictor variables µ interacting coefficient vector  Multilinear regression residuals xv Acknowledgements I first want to thank Dr. Sara Knox for the outstanding guidance and support you have given me from my first day as an undergraduate student researcher until now. You truly have made the largest impact on my academic development, and on this stage in my life, and I am forever grateful. I could not have been happier with my time in your lab. Of course, I couldn’t have done this without everyone at the Micrometorology Lab. Special thanks to Tin, for the guidance and company during my first foray into field work with the flux towers and data processing. To Nick, thanks for your great mentorship with the Manitoba sites and my research. To Zoran, thanks for teaching me so much about the flux systems and for all your patience required to do so. To Rick, thank you for showing me the ropes around the flux towers. I really value all the insights you’ve shown me in the field. To June, thanks for all the advice you’ve given me, whether it was with life or with academics, it was always helpful talking to you. To Zoe, thanks for showing me the other side of micromet field work, and for all the supportive talks in the office. And to everyone else in the lab that I’ve run out of space here for, thank you all for your help and camaraderie! Lastly, thank you to my committee members, Dr. Andy Black and Dr. Camilo Rey-Sanchez. Your enthusiasm and feedback for my research were beyond helpful. I’m especially grateful to Camilo, for so generously lending me your Flux Map code and all your help when I was just a hopeful undergrad. Additional thanks to the Natural Sciences and Engineering Research Council of Canada for funding parts of my undergraduate and Master’s research. xvi Dedication This thesis is dedicated to my incredibly supportive parents who had so much faith in me, my sister who is forever my role model, my girlfriend who kept me on my feet through this whole experience, and to my niece, for whom I hope to make a better world. 1 Chapter 1: Introduction The recent release of the Special Report on Global Warming of 1.5˚C (SR15) by the Intergovernmental Panel on Climate Change (IPCC) recommends a tightening of the Paris Climate Accords’ long-term temperature increase goal from 2˚C to 1.5˚C (IPCC, 2022). To achieve this goal with a limited overshoot of 0.5˚C, the second chapter of the SR15 states that we will have to halve our global carbon emission goal from 52-58 GtCO2eq yr-1 down to 25-30 GtCO2eq yr-1 (IPCC, 2022). This stricter emission target requires that we leverage as much of the natural climate change mitigating potential available to us as possible to help offset anthropogenic carbon emissions (Griscom et al., 2017; IPCC, 2022). Wetland ecosystems comprise 14% of total Natural Climate Solution opportunities around the globe, and play a significant role in the natural carbon cycle that removes carbon from the atmosphere (Griscom et al., 2017). In particular, wetlands are very effective at sequestering the carbon they remove from the atmosphere. By covering only about 5-8% of Earth’s land surface, they are capable of holding up to 30% of global soil carbon content, making them the densest terrestrial carbon stock among natural ecosystems (Kayranli et al., 2010; Mitsch & Gosselink, 2000; Petrokofsky et al., 2012). In Canada alone, an estimated 158 billion tonnes of carbon are stored in wetland soils, which is 60% more than the amount of carbon stored in Canadian forests (Roulet, 2000). This ecosystem service drives the growing need to advance our understanding of wetland carbon dynamics to ensure we can better utilize their carbon storage potential while avoiding mismanagement of their land-use that can lead to a net release of these dense carbon stocks back into the atmosphere (Bieniada & Strack, 2021; D′Acunha et al., 2019; Lee et al., 2017; McNicol et al., 2017). Methane gas (CH4) plays a crucial role in the wetland carbon cycle as it has a global warming potential 28-36 times that of CO2 sustained over 100 years, and its contribution to a 2 wetland’s carbon dynamics can determine whether the ecosystem ultimately has a net negative or positive climate forcing effect (Etminan et al., 2016; IPCC, 2022). 1.1 Wetland Carbon Cycle Atmospheric carbon is introduced into wetland ecosystems by photosynthesis from wetland vegetation (Kayranli et al., 2010; Roulet, 2000). Atmospheric carbon dioxide is drawn down and then converted into plant biomass and eventually soil organic carbon. Organic carbon can also enter (and leave) wetlands as either particulate organic carbon or dissolved organic carbon (Kayranli et al., 2010; Mitsch & Gosselink, 2000). The decomposition of organic matter via aerobic respiration, or heterotrophic respiration, into CO2 is the most energetically efficient use-path of terrestrial carbon, but is impeded by anaerobic conditions which arise from extended periods of high water levels (Griscom et al., 2017; Kayranli et al., 2010; Lee et al., 2017). Organic matter is able to accumulate in the soils when plant productivity exceeds the decomposition rate, leading to a net accumulation of organic matter in wetland soils (Mitsch & Gosselink, 2000). Draining wetlands, as was common 19th century practice in North America, introduces oxygen into the previously anoxic, carbon-rich wetland soils, and dramatically increases the release of carbon dioxide into the atmosphere through increased aerobic heterotrophic activity (Badiou et al., 2011; D′Acunha et al., 2019; McNicol et al., 2017; Mitsch & Gosselink, 2000). In undisturbed wetlands however, perennially high water levels create long-term anaerobic soils that support lower rates of heterotrophic respiration (D′Acunha et al., 2019; Kayranli et al., 2010; Mitsch & Gosselink, 2000; Torres-Alvarado et al., 2005). CO2 is also released via autotrophic respiration, and less energy-efficient ways of transforming the available carbon such as fermentation, and methanogenesis (Mitsch & Gosselink, 2000). 3 Autotrophic respiration is the biological conversion of carbohydrates into energy for cellular respiration in plants. This process uses atmospheric and soil O2 as its oxygen source. The frequent inundation of soils in wetlands however limits the availability of soil O2, and inhibits root respiration in the rhizosphere (Kayranli et al., 2010; Mitsch & Gosselink, 2000). The above-soil respiration creates CO2 as a by-product, and represents an important component of carbon efflux from wetlands (Mitsch & Gosselink, 2000). Fermentation is the process that breaks down high-molecular-weight carbohydrates into lower-molecular-weight acids and alcohols, and CO2 (Kayranli et al., 2010; Torres-Alvarado et al., 2005; Mitsch & Gosselink, 2000). The fermentative pathway taken is determined by the carbon matter available to be used as the terminal electron acceptor in the degradation process (Mitsch & Gosselink, 2000; Torres-Alvarado et al., 2005). For acids and alcohols, polysaccharides are used by fermentative bacteria to convert monomers into lactic acid and ethanol (Mitsch & Gosselink, 2000; Torres-Alvarado et al., 2005). The presence of other electron acceptors allow substrates to be broken down further into CO2 (Torres-Alvarado et al., 2005). The CO2 created by this process either diffuses back into the atmosphere through the soil and water column, used for anaerobic respiration, or used for methanogenesis as described in equation (1) (Mitsch & Gosselink, 2000). The most important outcome from this process is the creation of low-molecular-weight organic compounds which are central to anaerobic microbial activity (Mitsch & Gosselink, 2000; Wiebe et al., 1981). In wetlands, this is the primary source of carbon substrates used in methanogenesis (Torres-Alvarado et al., 2005). 4 1.2 Wetland Methane Emissions In methanogenesis, either CO2 or low-molecular-weight organic compounds are used as the terminal electron acceptor in the methanogen metabolism chemical reaction in equations (1) and (2) (Mitsch & Gosselink, 2000; Torres-Alvarado et al., 2005). 𝐶𝑂2 + 8𝑒− + 8𝐻+ → 𝐶𝐻4 + 2𝐻2𝑂 (1) 𝐶𝐻3𝐶𝑂𝑂𝐻 + 4𝐻2 → 2𝐶𝐻4 + 2𝐻2𝑂 (2) Equations (1) and (2) are just two examples of methanogen metabolism reactions that occur in wetlands. Boone et al. (1993) noted approximately 50 methanogen species within the Archaea domain, which can be grouped into four classes categorized by their carbon substrate metabolism (Torres-Alvarado et al., 2005). They all however share the common characteristics of not requiring oxygen for their reactions, and creating CH4 as a by-product (Torres-Alvarado et al., 2005). As a result, high methanogen activity and CH4 emissions are key characteristics that define the wetland carbon cycle (Bridgham et al., 2013). There are however other anoxic microbial reactions in wetlands that compete for carbon substrate and hydrogen ion resources. These reactions in descending order of redox potential are denitrification, manganese reduction, iron reduction, and sulfate reduction, with methanogenesis having the lowest potential for transformation among them (Mitsch & Gosselink, 2000; Torres-Alvarado et al., 2005). Because methanogenesis is the last stage of anaerobic degradation, the presence of nitrogen, manganese, iron, or sulfur in their oxidized forms (nitrate, manganic, ferric, and sulfate) will inhibit CH4 production, making methanogenesis and the accompanying gas emissions difficult to predict (Bridgham et al., 2013; Knox et al., 2021; Mitsch & Gosselink, 2000; Torres-Alvarado et al., 2005). Most freshwater wetlands, however, have low 5 concentrations of these higher-potential electron acceptors, making methanogenesis typically the predominating process in these wetlands (Torres-Alvarado et al., 2005). When unimpeded by its low redox potential, methanogen activity is primarily driven by the availability of simple carbon substrates, soil temperature (TS), and the volume of oxic vs. anoxic soil (Delwiche et al., 2021; Knox et al., 2021). Once created, CH4 is trapped in the soil until it is either oxidized and reduced, or transported and emitted into the atmosphere. CH4 oxidation occurs in the aerobic environments of the wetland by bacteria that use oxygen as the electron acceptor to convert CH4 to CO2 (Torres-Alvarado et al., 2005). This is a relatively slow process however, and requires the CH4 to interact with the aerobic environment for a long enough time for the bacteria to break down the gas (Torres-Alvarado et al., 2005). When production exceeds oxidation rates, CH4 that is not oxidized is emitted into the atmosphere via either diffusion, ebullition, or plant-mediated transport (Bridgham et al., 2013; Torres-Alvarado et al., 2005). Diffusion is the passage of CH4 along a concentration gradient in the soil, sediment, or water layer (Schipper & Reddy, 1994; Torres-Alvarado et al., 2005). This process is slow enough to allow some of the CH4 to oxidize before entering the atmosphere (Torres-Alvarado et al., 2005). Ebullition occurs when the CH4 concentration in the soil and sediment layer becomes large enough to overcome the hydrostatic pressure of the overlaying water layer to emerge through the sediment interface as gas bubbles (Bieniada & Strack, 2021; Torres-Alvarado et al., 2005). Plant-mediated transport occurs in aerenchymous vegetation that have adapted to allow gasses to ventilate through their stems and roots. Pressure gradients between the soil and plant tissue, created through either a build-up of gas concentration or convective flow, allow CH4 to diffuse into the plants, and then pass through the stem and stomata to enter the atmosphere (McNicol et 6 al., 2017; Torres-Alvarado et al., 2005). Ebullition and plant-mediated transport pathways typically produce much higher levels of CH4 emissions as they allow more of the gas to bypass oxidation processes (Bridgham et al., 2013; McNicol et al., 2017; Torres-Alvarado et al., 2005). The wetland carbon cycle is defined by the inflows and outflows of carbon dioxide and CH4 between land and atmosphere, and through lateral fluxes (Kayranli et al., 2010). In the conversation about wetlands as a nature-based climate solution, it is important to emphasize that CH4 is a greenhouse gas that has between 28-36 times the sustained global warming potential of CO2 over a 100-year time scale (Etminan et al., 2016; IPCC, 2022). Wetlands are also the single largest natural source of CH4 with global emission estimates of 115 to 227 Tg CH4 year-1, representing 25-40% of total natural and anthropogenic CH4 emissions (Bloom et al., 2010; Bridgham et al., 2013). Combining the strong global warming potential of CH4 with the high uncertainty arising from its spatio-temporal variability, it is critical to improve our understanding of wetland CH4 emissions to better quantify the role of wetlands as a climate mitigation solution. 1.3 Research objectives CH4 emissions are dependent on CH4 production, consumption, and release. Methanogenesis activity in the biosphere determines the CH4 pool in wetland soils, while oxidation and emission pathways controls how much and how quickly CH4 stocks are released into the atmosphere (Bridgham et al., 2013; Delwiche et al., 2021; Knox et al., 2021; Torres-Alvarado et al., 2005). Each of these processes are driven by a number of different biotic and abiotic factors, making the spatial distribution of CH4 emissions in wetlands extremely heterogeneous and difficult to predict (Bridgham et al., 2013; Knox et al., 2021; McNicol et al., 2017). This heterogeneity 7 makes it challenging to assess landscape-scale flux estimates from measurement methods like soil flux chambers and the eddy covariance (EC) method which only sample a portion of the broader ecosystem (Baldocchi, 2014; Chu et al., 2021; Forbrich et al., 2011; Oertel et al., 2016; Smith et al., 2003). Measurements sampled from areas of a wetland that are not representative of the ecosystem behaviour on the larger scale can lead to either under- or over-estimations of the ecosystem-scale fluxes (Chen et al., 2012; Chu et al., 2021). This is an especially difficult challenge with the EC method because the data is tied to a single flux tower’s location. These flux towers are expensive and can be very difficult to relocate, making it essential to have a good understanding of the area sampled by the flux tower and the representativeness of this area relative to the broader wetland ecosystem. A tower with sensors mounted at a measurement height of 2 metres above the vegetation surface will have an effective surface sampling area (flux footprint) of approximately 2 km2 as a general rule of thumb(Kljun et al., 2015; Lee et al., 2021). With some wetlands spanning up to 50 km2, a typical flux footprint area may only cover a small proportion of the whole ecosystem’s surface fluxes (D′Acunha et al., 2019). This leads to potential misrepresentation of CH4 surface features depending on how representative the flux tower footprint area is of the ecosystem CH4 features as a whole (Chen et al., 2012; Rey‐Sanchez et al., 2022). It is therefore important to advance our understanding of the spatial distribution of wetland methane fluxes (FCH4), and to develop approaches to better analyze and interpret the spatial heterogeneity and their drivers. Previous studies have either combined eddy covariance and remote sensing to investigate the spatial representativeness of flux towers from a land cover change or vegetation characteristics perspective, or have investigated the spatial behaviour of CH4 emissions using chamber 8 measurements or multi-tower approaches (Bieniada & Strack, 2021; Chen et al., 2012; Chu et al., 2021; D′Acunha et al., 2019; Lu et al., 2017; McNicol et al., 2017). However, no study has yet investigated the within-footprint spatial variability of CH4 emissions from individual eddy covariance flux towers. This study explores a novel method of characterizing the within-footprint spatial variability of CH4 emissions in wetlands using a combination of Footprint-weighted Flux Maps and remote sensing. This is a top-down approach to analyzing the surface conditions within an eddy covariance flux tower’s flux footprint, and identifying the sources and environmental drivers of methane heterogeneity within the footprint. The most significant of these drivers as detectable by remote sensing are: • Gross primary productivity (GPP) – Biological control of substrates required for methanogenesis, and biological control for plant-mediated methane transport. • Soil temperature (TS) – Physical control of environmental conditions for microbial processes. • Water table depth (WTD) – Physical control of anaerobic soil zone thickness, and physical control of soil saturation and compression which determines methane diffusion and ebullition. (Bridgham et al., 2013; Delwiche et al., 2021; Knox et al., 2021; Smith et al., 2003) This spatial characterization method will be used to address two objectives in this thesis: 1) Determine whether there exist significant relationships between wetland surface conditions within a flux footprint and FCH4. This is explored by investigating 9 heterogeneous FCH4 features using Footprint-Weighted Flux Maps, identifying distinct surface conditions through Landsat 8 remote sensing indices, and then testing the correlation relationships and strengths with simple linear regression. 2) If statistically significant correlations exist, extrapolate and upscale the relationships to estimate FCH4 beyond the reach of EC flux towers. This is explored with multiple linear regression modelling. The scope of this study is limited to freshwater wetlands as these ecosystems exhibit stronger CH4 emissions than saltwater wetlands, which typically have very low CH4 emissions due to the inhibition of methanogenesis by sulfates present in saltwater (Lu et al., 2017; Mitsch & Gosselink, 2000). 10 Chapter 2: Methods 2.1 Study sites Four freshwater wetland sites were selected for this study: Hogg Marsh (MBPPW1), Young Marsh (MBPPW2), Mayberry Marsh (US-Myb), and Winous Point North Marsh (US-WPT). Eddy covariance data from Hogg and Young were collected by Ducks Unlimited Canada in partnership with The University of British Columbia, while data from US-Myb and US-WPT were downloaded from the Ameriflux database. These sites, with the exception of MPBBW1, were chosen for having high spatial heterogeneity of land surface conditions. The sites also had to have records of Monin-Obukhov length, friction velocity, wind direction, and wind speed variances, which are essential atmospheric turbulence inputs for the Footprint-weighted Flux Maps. Their data also had to be available during Landsat 8’s operational period from February 2013 to present. MBPPW1, located at 50.3705˚N, -100.5339˚W, is an isolated grassland marsh in the Prairie Pothole Region of Canada (Figure 1). The land cover at this site is fairly homogenous and dominated by emergent vegetation, primarily populated by Schoenoplectus tabernaemontani (scirpus) and Typha (cattails) spp. The water chemistry in this wetland is characterized by high sulphate concentrations, which create conditions of low CH4 emissions through the impediment of methanogenesis (Soued et al., in preparation). FCH4 at this site is therefore very low with only slight increases in CH4 emissions in the growing seasons, as seen in the time series plot in Figure 5a. 11 Figure 1 MBPPW1. Photograph on left is a satellite image of MBPPW1 taken from Google Earth. The yellow star marks the flux tower location. Photograph on right is a photo of the site taken on 2022/07/18. MBPPW2, at 50.3623˚N, -100.20242˚W coordinates, is an isolated cropland marsh also located in the Prairie Pothole Region of Canada (Figure 2). Unlike MBPPW1, the land surface at this site is extremely heterogeneous with a combination of open water and emergent vegetation. Large mats of submersed macrophytes are found near the open water surface during the growing seasons, and the emergent vegetation is predominantly populated by Typha spp. Low water salinity at this site, enables high CH4 production and strong seasonal CH4 emissions in the warmer months (Figure 5b). 12 Figure 2 MBPPW2 Photograph on left is a satellite image of MBPPW2 taken from Google Earth. The yellow star marks the flux tower location. Photograph on right is a photo of the site taken on 2021/05/19. US-Myb is a restored, emergent freshwater marsh located in the western region of the Sacramento-San Joaquin Delta in California, United States (Schipper & Reddy, 1994). The hydrology at this site is influenced by rainfall and water pumped in from an adjacent river. The land surface surrounding its flux tower is characterized by a mix of open water with macrophyte mats and emergent vegetation populated by Schoenoplectus spp. and Typha spp. Figure 3 US-Myb Photograph on left is a satellite image of US-Myb taken from Google Earth. The yellow star marks the flux tower location. Photograph on right is a photo of the site taken on 2011/10/07. 13 US-WPT is an isolated marsh in northern Ohio, United States characterized by open water and emergent vegetation. The open water vegetation is dominated by Nymphaea odorata (water-lilies) and Nelumbo lutea (yellow lotus) while the emergent vegetation is populated by Typha angustifolia, Hibiscus moscheutos, and Sparganium americanum (Chu et al., 2016). Figure 4 US-WPT Photograph on left is a satellite image of US-WPT taken from Google Earth. The yellow star marks the flux tower location. Photograph on right is a photo of the site taken on 2011/06/16. 14 Figure 5 Methane flux time series. Time series plots of half hourly methane fluxes (nmol/m2/s) from each of the four study sites. Subplot a) is the time series for MBPPW1, b) for MBPPW2, c) for US-Myb, and d) for US-WPT. Each subplot has individualized y-ranges to capture the varying CH4 flux ranges at each site. 2.2 Eddy covariance data Half-hourly gas, momentum, and energy fluxes between the surface and the atmosphere were measured using eddy covariance (EC) sensors mounted on in-situ wetland flux towers. Cleaned, 15 high quality EC data were downloaded from the Ameriflux database for US-Myb and US-WPT, while EC data for MBPPW1 and MBPPW2 were collected by the Institute of Wetland and Waterfowl Research at Ducks Unlimited Canada, and cleaned and gap-filled by the Micrometeorology Lab at the University of British Columbia. 2.2.1 Eddy covariance measurements at MBPPW1 and MBPPW2 At MBPPW1 and MBPPW2, high-frequency 20 Hz data were measured with an open-path infrared gas analyzer (IRGA) (Model LI-7500DS, LI-COR Inc., Lincoln, NE, USA) for CO2 and H2O molar densities, an open-path IRGA (Model LI-7700, LI-COR Inc.) for CH4 molar density, and a 3-dimensional ultrasonic anemometer (Model 81000, R.M. Young Co., Traverse City, MI, USA) for wind speed and wind direction in 3 dimensions (u, v, and w). CO2, H2O, and CH4 fluxes were then calculated by Reynold’s averaging and the product of vertical wind fluctuations and trace gas’s molar density (Burba & Anderson, 2010). These sensors were mounted on flux towers at measurement heights of 4.08m for MBPPW1 and 4.09m for MBPPW2. Data were logged using LI-COR’s Synchronization, Management and Real Time Flux system (Smartflux 3) software and onsite data pre-processing was done using a LI-COR 7550 Analyzer Interface Unit. Final re-processing of the high-frequency data was done using the EddyPro v7.0.6 software. De-trended with half-hourly block-averaging and corrections for gas concentration sensitivity to air temperature and water vapour using Webb-Pearman-Leuning equations were performed before passing the data to post-processing steps to filter and gap-fill fluxes (Burba & Anderson, 2010. 16 2.2.1.1 Data filtering Filtering for the Manitoba sites was done in four stages. The first stage removed data that fell outside of realistic minimum and maximum ranges. The ranges specified for each meteorological variable were informed either by physical limitations, such is the case for relative humidity, or by regional climate, such as for variables like air temperature (TA), incoming shortwave radiation (SWIN), and wind speed. Realistic ranges for more site-specific variables such as sensible heat, latent heat, CO2, H2O, and CH4 fluxes were determined by observation. The second stage of filtering removed data coming from poor sensor conditions. These are conditions when Smartflux system voltage fell below 0 V or rose above 30 V, Data Retention Module voltage fell below 0V or above 50V, Data Acquisition Module temperature fell outside of -60 to 60˚C, and when LI7500 and LI7700 signal strengths fell below 20%. An additional pass of filtering removed flux data from measurements when wind directions are between 330˚ and 30˚. This was done to remove turbulence interference from the flux towers which are positioned north of the EC sensors. The third stage diagnosed half-hourly flux data quality by stationarity and atmospheric turbulence with quality tests proposed by Foken at al. (2004). In this stage, the poorest quality data were given a “2” flag, data that were suitable for general analysis were given a flag of “1”, and the highest quality data were given a flag of “0”. The methods in this study benefitted from larger sample sizes while also requiring an adequate standard of quality in the data, so the Foken flag threshold for this study was set at 1, and all data with a flag of “2” were removed. 17 The fourth and final stage of data filtering removed bad data that resulted from insufficient atmospheric turbulence. Because the EC method relies on the covariance between gas concentrations and vertical wind speeds to measures gas fluxes (Baldocchi, 2014), adequate atmospheric turbulence is required to produce high quality flux calculations. This turbulence can be represented by friction velocity, a variable that measures vertical shear stress in units of metres per second. When conditions are too stable and friction velocity is low, the EC method will tend to underestimate fluxes and misrepresent real micrometeorological processes (Barr et al., 2013). A minimum friction velocity threshold was therefore applied to remove EC measurements under insufficient turbulence conditions. Friction velocity thresholds were estimated with moving point tests via the REddyProc v 2.2.0 R package (Wuzler et al., 2018). Half-hourly flux measurements for which friction velocities fell beneath the threshold were removed from the final dataset. 2.2.1.2 Data gap-filling Three different approaches to gap-filling were employed for three separate variable classes in the dataset: meteorological variables, CO2 and heat fluxes, and CH4 fluxes. 2.2.1.2a Gap-filling meteorological data Meteorological conditions are expected to influence large geographical areas so there is high confidence that neighbouring climate stations should have very similar half-hourly observations. If measurements have high similarity, then gaps in the dataset can be filled with climate station data. 1:1 Pearson correlation tests on major meteorological variables were performed between MBPPW1 and MBPPW2 and their respective nearby climate stations to check the similarity in 18 observations. MBPPW1 had two nearby climate stations, called Hamiota and Oakburn, so the correlation test was also used to choose which station was more appropriate to fill data gaps from (see Appendix A for detailed methods). MBPPW1 was gap-filled using Oakburn climate station located 26 km north of the site, and MBPPW2 was gap-filled using Newdale climate station located 4 km north west of the site. 2.2.1.2b Gap-filling CO2 and heat flux data CO2, sensible heat (H) and latent heat (LE) fluxes, which are strongly controlled by site-level environmental conditions, were gap-filled with the Marginal Distribution Sampling (MDS) technique using the REddyProc v 2.2.0 package in R (Reichstein et al., 2005; Wutzler et al., 2018). Missing measurements in these fluxes were modelled with data calculated by their correlation with their driver variables: Net radiation, TA, and vapour pressure deficit. 2.2.1.2c Gap filling FCH4 FCH4 are much more complicated to model due to their non-stationary behaviour and temporally non-uniform relationships with their driver variables (Knox et al., 2021; Kim et al., 2002). Random forest modelling was used to gap-fill CH4 following the methods described in Kim et al. (2002). The predictor variables used to train the random forest models were friction velocity, net ecosystem exchange (NEE), LE, H, incoming shortwave radiative flux, TA, relative humidity, vapour pressure deficit, day of year, wind direction, and TS at 40 cm above, 10 cm above and 10 cm below the sediment-water interface. The random forest variable importance graph for MBPPW1 (Figure 6) shows that H, NEE, and TS 10 cm below sediment-water interface were found to be the strongest drivers of FCH4 at MBPPW1. 19 Figure 6 (MBPPW1) Random forest driver importance graph. Methane driver variable inputs are listed on the y-axis, and their relative importance to methane prediction are plotted on the x-axis. TS 10 cm below the sediment-water interface, wind direction, day of year, and H flux were the most dominant driver variables in the model for MBPPW2 (Figure 7). 20 Figure 7 (MBPPW2) Random forest driver importance graph. Methane driver variable inputs are listed on the y-axis, and their relative importance to methane prediction are plotted on the x-axis. 2.3 Footprint-weighted Flux Maps Footprint-weighted Flux Maps following the flux heterogeneity mapping technique described in Rey-Sanchez et al. (2022) were used to find the spatial distribution of CH4 fluxes within each site’s flux footprint climatology. This approach uses half-hourly atmospheric turbulence data over a specified period of time to generate flux footprint iterations, and superimposes them on a surface grid to create a footprint climatology. A flux footprint is a transfer function that quantifies the importance of a certain position in a land surface to the measurement made at the flux tower (Kljun et al., 2015; Kormann & Meixner, 2001). The eddy covariance flux density 21 measured at height zm is calculated by Equation 3, where FEC(0,0,zm) is the vertical flux of the scalar of interest, F(x,y,0) is the surface flux at x and y, and f(z,y,zm) is the footprint function (Kormann & Meixner, 2001; Rey‐Sanchez et al., 2022). 𝐹𝐸𝐶(0,0, 𝑧𝑚) = ∫ ∫ 𝐹(𝑥, 𝑦, 0)𝑓(𝑥, 𝑦, 𝑧𝑚)𝑑𝑥𝑑𝑦∞0∞−∞ (3) The footprint function f(z,y,zm) is calculated by the product of the cross-integrated footprint function fy and the crosswind Guassian plume dispersion function Dy (Equation 4) 𝑓(𝑥, 𝑦, 𝑧𝑚) = 𝑓𝑦(𝑥, 𝑧𝑚)𝐷𝑦(𝑥, 𝑦) (4) 𝐷𝑦(𝑥, 𝑦) =exp (−𝑦22𝜎𝑦2)√2𝜋𝜎𝑦 (5) where y is the standard deviation of the plume in the y dimension along the stream-wise x direction (Amiro, 1998; Kormann & Meixner, 2001; Rey‐Sanchez et al., 2022). Half-hourly FCH4 are then weighted by the footprint-grid intersections to calculate the average flux measured at a given grid cell (Rey-Sanchez et al., 2022). This weighting is calculated by: 𝐹𝑤𝑖,𝑗 =∑ 𝐹𝐸𝐶𝑚𝑓𝑖,𝑗,𝑚𝑀𝑚∑ 𝑓𝑖,𝑗,𝑚𝑀𝑚 (6) where M is the number of half-hourly points in a given observation period, FECm is the EC flux at time m, and 𝑓𝑖,𝑗,𝑚is the value of the footprint function at cell i,j, at time m (Rey‐Sanchez et al., 2022). The resulting grid provides a land-cover-independent spatial map of relative fluxes around a flux tower that is much more capable of detecting hotspots of strong CH4 emissions than previous footprint-aggregation approaches such as those described by Amiro (1998) and Chen et al. (2009). The flux footprint model used in this study was the Korman and Meixner (2021) analytical model (Equation 4). This footprint model was chosen over other commonly used models such as Kljun (2015) and Hsieh (2000), as it was found from tracer experiments to perform better in the detection of trace gas hotspots under highly heterogeneous conditions (Rey‐ 22 Sanchez et al., 2022). The footprint contours were limited to 80% to preserve confidence in the maps’ signal-to-noise ratio at the spatial level without over-sacrificing total footprint area coverage (Chu et al., 2021; Rey‐Sanchez et al., 2022). One Flux Map for each observation in this study was generated from a four-week period, centred on a Landsat 8 image collection date. The Flux Map datasets represented the FCH4 features surrounding each flux tower, and the corresponding land-surface conditions are represented by the Landsat 8 remote sensing datasets. 2.4 Landsat 8 Remote Sensing Remote sensing data in this study were taken from the United States Geological Survey Landsat 8 satellite. Four remote sensing products were chosen as proxy measurements of variables that are known to influence CH4. These products include: NDVI, NDWI, NDMI, and Land Surface Temperature (LST). They represent proxies for GPP, water table height, moisture conditions, and TS, which are all known controls of FCH4 (Bridgham et al., 2013; Delwiche et al., 2021; Knox et al., 2021; Oertel et al., 2016; Schipper & Reddy, 1994). (Table 1). 23 Table 1 Summary table of CH4 flux drivers and the associated remote sensing products used to approximate them. CH4 flux driver Reference Remote sensing product Wavelength band calculation Reference GPP (Colmer, 2003; Knox et al., 2021; Turetsky et al., 2014) NDVI 𝜆𝑁𝐼𝑅 − 𝜆𝑅𝑒𝑑𝜆𝑁𝐼𝑅 + 𝜆𝑅𝑒𝑑 Chen et al., 2012; Sellers, 1985 Water table height (Knox et al., 2021; Oertel et al., 2016; Schipper & Reddy, 1994) NDWI 𝜆𝐺𝑟𝑒𝑒𝑛 − 𝜆𝑁𝐼𝑅𝜆𝐺𝑟𝑒𝑒𝑛 + 𝜆𝑁𝐼𝑅 Rokni et al., 2014; McFeeters et al., 1996 Moisture conditions (Bridgham et al., 2013; Knox et al., 2021; Schipper & Reddy, 1994) NDMI 𝜆𝑁𝐼𝑅 − 𝜆𝑆𝑊𝐼𝑅𝜆𝑁𝐼𝑅 + 𝜆𝑆𝑊𝐼𝑅 Rokni et al., 2014; Gao, 1996 Soil temperature (Knox et al., 2021; Torres-Alvarado et al., 2005) LST 𝐴𝑖𝑇𝑏𝜀+ 𝐵𝑖1𝜀+ 𝐶𝑖 Ermida et al., 2020 Water table height drives methanogenesis through its control on anoxic soil volume and methanogen activity (Knox et al., 2021; Schipper & Reddy, 1994). Water table height can be measured either as the degree of open water or water saturation in the rhizosphere (Gao, 1996; Rokni et al., 2014). Open water can be approximated by the Normalized Difference Water Index (NDWI), calculated as the difference in apparent reflectance between the green band (0.53-0.59 µm) and the near infrared band (0.85-0.88 µm) (Rokni et al., 2014; McFeeters et al., 1996) divided by their sum. Soil and vegetation moisture content can be approximated by the Normalized Difference Moisture Index (NDMI), calculated as the difference between the apparent reflectance in the near infrared band and the shortwave infrared band (1.57-1.65 µm) divided by their sum (𝜆𝑁𝐼𝑅−𝜆𝑆𝑊𝐼𝑅𝜆𝑁𝐼𝑅+𝜆𝑆𝑊𝐼𝑅) (Gao, 1996; Rokni et al., 2014). 24 Gross primary production (GPP) is a biological control for both CH4 production as well as CH4 release. Root exudation of organic carbon into the rhizosphere from productive vegetation provides fuel for methanogens (Knox et al., 2021; Turetsky et al., 2014). Dominant wetland vegetation adapt to inundated conditions by developing aerenchymous tissue to deliver oxygen to their roots, but this also creates a pathway for CH4 to pass from the soil into the atmosphere (Colmer, 2003). Vegetation greenness is a strong indicator of photosynthetic capacity and GPP, and can be quantified with the Normalized Difference Vegetation Index (NDVI), calculated by the difference in apparent reflectance between the near infrared and red wavelength band (0.64-0.67 µm) divided by their sum (Chen et al., 2012; Rokni et al., 2014; Sellers, 1985). TS is a physical control for methanogen activity. There is a direct relationship between temperature and microbial activity (Torres-Alvarado et al., 2005). This, along with GPP, are the primary drivers for seasonality in CH4 emissions in temperate wetlands, typically stimulating higher CH4 fluxes in the growing season (Knox et al., 2021; Torres-Alvarado et al., 2005). Although, it should also be noted that increased temperatures also increase the CH4 consumption activity of methanotrophs (Torres-Alvarado et al., 2005). TS is approximated by Land Surface Temperature (LST), which is calculated in this study using the open-source Google Earth Engine code from Ermida et al. (2020). This code utilizes the Top of Atmosphere (TOA) brightness temperatures from Landsat 8’s thermal infrared wavelength bands (10.6-11.19 µm), resampled to 30 m spatial resolution (Ermida et al., 2020). LST is calculated by 𝐴𝑖𝑇𝑏𝜀𝑠+ 𝐵𝑖1𝜀𝑠+ 𝐶𝑖, where the surface emissivity (𝜀𝑠) at each pixel is calculated by the Vegetation-Cover Method, derived from fractional vegetation cover and fractional difference between bare and vegetated soil NDVI 25 values, Tb is the TOA brightness temperature, and Ai, Bi, and Ci are radiative transfer coefficients (Ermida et al., 2020). All remote sensing data were collected from the Landsat 8 tier 1 collection via Google Earth Engine to use the highest quality data available from the satellite (Gorelick et al., 2017). Each observation consists of four Landsat 8 images, visualizing each of the study’s remote sensing products at a 30-m resolution, and centred on each site’s EC flux tower. Available Landsat 8 images were limited by the satellite’s 16-day orbital period, as well as cloud and snow coverage over the flux tower site. Observations in this study were taken from periods that had unobstructed Landsat 8 images available during the growing season when CH4 spatial variability would be the most pronounced, and that had fairly even samplings of wind direction to help lower biases in land-area coverage. 2.5 Map data processing Each observation comprised a Footprint-weighted Flux Map created from four weeks of EC data, and a Landsat 8 image consisting of 4 observation variables: NDVI, NDWI, NDMI, and LST. Before these two datasets were studied together, they were processed to have matching spatial scales. This entailed isolating Landsat 8 maps’ spatial coverage to only the footprint area within by the Flux Map climatologies, and then matching the spatial resolution of the finer-scale Flux Maps to the coarser-resolution Landsat 8 maps. Finally, to be able to better isolate and compare surface features between the two maps, hierarchical clustering was applied to both datasets. 26 2.5.1 Matching spatial layout and coverage Landsat 8’s spatial coordinates in latitude and longitude were first converted into distance from the flux tower in metres, where the flux tower was located at the (0, 0) coordinates. The edge limits of the new Landsat maps were defined by the boundaries of the Flux Map footprint climatologies at the 80% footprint contour, which were also measured in metres from the tower. The Landsat 8 pixels, while accounting for their 30-metre resolution, were trimmed to remove any data points with at least 50% of their area outside of the footprint edge limits. An example of a Landsat 8 map matched to its corresponding Flux Map’s footprint area for MBPPW2 is shown in Figure 8. Figure 8 Footprint-matched Landsat 8 map. An example Landsat 8 map of NDVI matched to a Flux Map footprint area. Figure 8a shows the raw NDVI image collected from Google Earth Engine. Figure 8b shows the Landsat 8 map after removing pixels outside of the footprint area. The black lines mark the horizontal and vertical coordinates of the flux tower. 2.5.2 Matching pixel resolution The Footprint-weighted Flux Maps have a 4-m spatial resolution, whereas the Landsat 8 pixels have a 30-metre resolution. To match the two properties, Flux Map resolutions was coarsened to 30 metres by grouping the Flux Map cells into 30 by 30 metre bins centred on the coordinates from each Landsat 8 pixel. The bins were then represented by their average flux values. An example taken from MPBBW2 is shown in Figure 9. The final Flux Map datasets then had the 27 same number of data points as the Landsat 8 maps, and each flux data point had a pixel counterpart in the Landsat maps that occupied the same spatial coordinates. With matching map areas and spatial resolutions, the Flux Map and Landsat 8 datasets could then now be analyzed together, but before proceeding to regression analyses, one last pre-processing step was needed to address spatial artifacts found in the Flux Maps that can introduce potential spatial errors to the flux datasets. This pre-processing step was performed with hierarchical clustering. Figure 9 Example Flux Map resolution coarsening. An example of coarsening a Flux Map resolution to match Landsat 8’s 30 m resolution taken from MBPPW2. a) is the original Flux Map at 4-m resolution, and b) is the 30-m-binned Flux Map. 2.5.3 Hierarchical clustering Spatial artifacts are created when features of pronounced FCH4 behaviour are located at some distance from the flux tower at origin (0, 0) and are detected by flux footprints, but the radial nature of footprint calculations conflates the feature FCH4 with measurements along the range between the feature location and (0, 0). These artifacts are seen as residual “shadows” of the flux feature’s characteristics cast back to the origin point. These shadows obscure the underlying pixels’ actual FCH4 information (Figure 9a). This causes what should be isolated pockets of hotspots or low emissions to appear as long slices of high/low fluxes stretching from the origin to the feature location (Rey‐Sanchez et al., 2022). To reduce their impact on the datasets, hierarchical clustering was used to group Flux Map data points into clusters of distinct flux 28 features. A bottom up, agglomerative clustering approach with Ward’s linkage method was used to create these clusters by identifying data points, that when added to a cluster group, minimize the increase in total within-cluster variance (Gong & Richman, 1995). This clustering approach was chosen over other methods for this study as it was found to perform more accurately with relatively small sample size geospatial data (Gong & Richman, 1995). To determine the optimal number of cluster groups, dendrogram analyses were performed to identify the cut-off points (Figure 10). Finding an appropriate cut-off point is important. Having a large enough number of cluster groups is crucial to allowing all distinctly significant map features to be captured, while a number that is too large risks breaking up features into non-distinct noise. 29 Figure 10 Hierarchical clustering dendrogram analysis. Figure 10a shows a Footprint-weighted Flux Map generated from US-WPT from June 5 to July 5, 2023. Figure 10b) shows the dendrogram generated from Ward’s linkage hierarchical clustering. Figure 10c) shows the clustered Flux Map with data points colour-coded to indicate their cluster group. The dendrograms show the hierarchical organization of data points grouped down the graph into branches of increasing cluster numbers. Each limb of the dendrogram represents one cluster group. The cut-off point is determined by reading the graph from top to down, and finding the 30 last branch, beneath which every subsequent branch is tightly grouped and separated by only small distances. The cut-off cluster number is then the number of limbs intersected by a horizontal line drawn at this last branch. In the sample dendrogram shown in Figure 10b, this cut-off branch is the first of the orange-coloured branches, and the horizontal red line denotes the cut-off point. Five limbs intersect the cut-off line, indicating the optimal number in this example to be five cluster groups. This assessment can be confirmed with a visual evaluation of the Footprint-weighted Flux Map in Figure 10a and the clustered Flux Map in Figure 10c to confirm there are indeed five distinct flux features. Based on multiple tests, the optimal cluster number with these Flux Maps was discovered to be at a minimum of four, with most observations being able to accommodate five clusters, and only few instances of Flux Maps being able to fit up to six clusters. Five clusters were then set as the default parameter for all observations in this study with the aim reducing noise while maximizing available data points. Dendrogram analyses were always performed with each observation to ensure there was no over-fitting. After the Flux Maps were clustered, the same cluster groupings were applied to the Landsat 8 images. Each map then had five clusters of remote sensing data, spatially aggregated in the same way as the Flux Maps. The average cluster values were then calculated for each cluster, and then collected to create a final dataset comprising five data points for each of the maps. 2.6 Linear regression analysis Simple linear regression paired with Pearson correlation tests were performed on the clustered Landsat 8 and Footprint-weighted Flux Map datasets to determine the strength and significance of their relationships. Individual linear regression tests were performed between each of the four 31 remote sensing products and clustered FCH4. To test for significance, Pearson correlation tests, using the the scipy 1.10.1 pearsonr python package, were taken at a significance level of α = 0.10 (Zwillinger & Kokoska, 2000). Pearson correlation coefficients (r-value), ranging from -1 to 1, determined the strength and direction of the relationships between FCH4 and each remote sensing product. Regressions with Pearson p-values larger than the 10% significance level were excluded from the final analyses. 2.7 Multiple linear regression modelling Multiple linear regression (MLR) models, fit using the lm function in the stats R package version 3.6.2 (Chambers et al., 1973), were used to capture the specific relationships between FCH4 and the four Landsat 8 predictor indices at MBPPW2. Here we decided to focus on one site, MBPPW2, to illustrate the potential of leveraging remote sensing products to upscale FCH4 in heterogeneous wetlands. NDVI, NDWI, NDMI, LST, and month-of-the-year (month.f) were used as fixed effect predictor variables, and FCH4 was the dependent variable. Month.f was included to account for temporal trends in predictor and dependent variables over the growing season. Clustered and footprint-matched Landsat 8 index data and Flux Map FCH4 data from MBPPW2 over 2021 to 2022 were used to train the MLR models. Data from two observations, 2021-07-14 and 2021-09-25, were excluded from the training dataset, and used for testing. Two MLR models were tested in this study. The first was a model with all the available variables (NDVI, NDWI, NDMI, LST, and month.f) and no interaction terms between them (Simple Model). The second MLR model was found using an in/out stepwise selection algorithm starting 32 with all predictor variables and their interaction terms (Model Stepwise). Model selection was based on achieving the highest adjusted R2 metric. The two MLR models can be expressed as: 𝐹𝐶𝐻4 = 𝑋𝛽 + 𝑍𝜇 + 𝜀 (7) 𝜀 ~ 𝑁(0, 𝜎𝜀2) (8) where X is the matrix of respective predictor variables, β is the predictor coefficient vector, Z is the matrix of interacting predictor variables with month, and µ is the interacting coefficient vector, and ε are the independent and identically distributed residual terms. The Shapiro-Wilk Normality Test was used to ensure that the residual normality assumption was met for both models. Linearity and residual homoscedasticity assumptions were assessed visually with predicted-y, and residual plots respectively. Model results from the training data were tested for significance at the 10% significance level to be consistent with the linear regression methods. 33 Chapter 3: Results The simple linear regression analyses were able to find strong correlations with each of the Landsat 8 remote sensing indices. Among the statistically significant observations across all four sites, mean Pearson r in absolute values found for each index was 0.916 for NDVI, 0.892 for NDWI, 0.898 for NDMI, and 0.945 for LST. The high correlation strengths suggest a strong potential for these methods to be applied in FCH4 spatial-variability characterization, and FCH4 upscaling and prediction. 3.1.1 Simple linear regression – MBPPW1 Relationships between FCH4 and each of the clustered Landsat 8 products (NDVI, NDWI, NDMI, and temp) are illustrated in scatter plots (2021 shown in Figure 11, and 2022 shown in Appendix Figure 4). Each observation is identified by its Landsat 8 image acquisition date, and the linear regressions for its Landsat 8 indices are grouped by row. Relationships that did not meet the 10% significance level are highlighted in red, and excluded from the final analysis. Of the 40 linear regression relationships calculated at MBPPW1 (four indices over 10 observations), 16 of them were statistically significant at the 10% significance level. These relationships demonstrated strong correlation strengths with the average Pearson r in absolute values for NDVI, NDWI, NDMI, and LST being 0.862, 0.855, 0.898, and 0.962 respectively. However, the signs of these relationships were not consistent across observations. Interpreting negative r-values as negatively correlated relationships between two variables, and positive r-values as positively correlated relationships, three of the regressions for NDVI were negative on May 2021, July 2022, and August 2022, while two were positive on July and August 2021 34 (Table 2). Similarly, for NDWI, two of the five significant observations had negative correlations, and the other three were positively correlated. For NDMI and LST, significant observations were split into equal halves of negative and positive relationships. A boxplot showing the distribution of r-values and their medians (Figure 12) show that across 2021 and 2022, FCH4 tended to relate negatively with NDVI and NDMI, while it tended to relate positively with NDWI. 35 Figure 11 MBPPW1 simple linear regression. Scatterplots of simple linear regression analyses for MBPPW1 over the 2021 period. Each row is an individual observation, identified by the Landsat 8 image acquisition date on the left. Highlighted in red are Pearson correlation values of runs whose p-values are larger than 0.1. 36 Table 2 MBPPW1 r-value data table. Statistically significant Pearson correlation coefficient values for MBPPW1 assessed at the α = 0.1 significance level. Highlighted in orange are negative coefficient values. The double line separates the results from 2021 and 2022. Observations where all correlations were not significant are omitted. MBPPW1 Observation NDVI NDWI NDMI LST 2021-05-11 -0.938 0.832 0.927 - 2021-05-27 - -0.886 -0.936 - 2021-06-12 - - - 0.962 2021-07-14 0.881 -0.876 - - 2021-08-31 0.817 - - - 2022-07-01 -0.834 0.837 -0.844 - 2022-08-02 -0.830 0.847 0.884 - 2022-09-03 - - - -0.963 Figure 12 Boxplots of statistically significant Pearson correlation coefficients for MBPPW1 Data are separated into 2021 and 2022 periods in the left two subplots, and the right-most plot is created from the full 2021 - 2022 dataset. Correlation coefficient data points are plotted as blue circles, and medians are plotted as the red lines. 3.1.2 Simple linear regression – MBPPW2 Simple linear regression scatterplots for 2021 (Figure 13) show that 13 observations were made at MBPPW2 with 22 statistically significant regressions found, and for 2022 (Appendix Figure 5), 9 statistically significant regressions were found. The mean Pearson r in absolute values was 37 0.931, 0.899, 0.892, and 0.975 for NDVI, NDWI, NDMI, and LST respectively. Compared to MBPPW1, the relationship signs at this site were noticeably more consistent across observations. From the r-values collected in Table 3, NDVI’s relationship with FCH4 was consistently negative across 2021 and 2022 with the exception of one observation on 2022-07-26. NDWI and NDMI were both consistently positively related to FCH4. LST was slightly more variable with one observation showing a positive relationship in July 2021, and two observations showing negative relationships on 2022. The median r-values, shown in a boxplot (Figure 14) across the two years was -0.925 for NDVI, 0.922 for NDWI, 0.865 for NDMI, and -0.897 for LST. 38 Figure 13 MBPPW2 simple linear regression. Scatterplots of simple linear regression analyses for MBPPW2 over the 2021 period. Each row is an individual observation, labelled by the Landsat 8 image acquisition date on the left. Highlighted in red are Pearson correlation values of runs whose p-values are larger than 0.1. 39 Table 3 MBPPW2 r-value data table. Statistically significant Pearson correlation coefficient values for MBPPW2 assessed at the α = 0.1 significance level. Highlighted in orange are negative coefficient values. The double line separates the results from 2021 and 2022. Observations where all correlations were not significant are omitted. MBPPW2 Observation NDVI NDWI NDMI LST 2021-07-14 -0.977 0.972 0.975 0.975 2021-08-08 -0.925 0.899 0.865 - 2021-08-31 - 0.809 - - 2021-09-09 -0.891 0.867 0.836 - 2021-09-25 -0.931 0.950 - - 2022-05-23 -0.95 0.922 - -0.991 2022-07-26 0.883 - - - 2022-09-03 -0.899 0.961 0.980 -0.897 2022-09-28 - - 0.809 - Figure 14 Boxplots of significant Pearson r-values for MBPPW2. Data are separated into 2021 and 2022 periods in the left two subplots, and the right-most plot is created from the full 2021 - 2022 dataset. Correlation coefficient data points are plotted as blue circles, and medians are plotted as the red lines. 3.1.3 Simple linear regression – US-Myb Results from simple linear regression analyses on 2021 data from US-Myb are plotted in Figure 15. Thirteen of the 32 regressions were found to be significant at the 10% significance level. Eight of the insignificant regressions failed the significance test by only a slim margin with p-values falling between 0.1 and 0.15. Similar to the other two sites, the significant regressions at 40 US-Myb exhibited strong correlation coefficients with mean Pearson r in absolute values of 0.942, 0.877, 0.908, and 0.921 for NDVI, NDWI, NDMI, and LST respectively. The significant r-values listed in Table 4, and their distribution is shown in the boxplot in Figure 16. FCH4 exhibited consistently positive relationships with NDMI, and consistently negative relationships with LST. The removal of high p-value measurements left equal numbers of regressions on either side of the r = 0 line for NDVI, and NDWI (Figure 16). No clear trends can be determined for NDVI and NDWI at the α = 0.1 significance level. At the α = 0.15 significance level however (summarized in Appendix Table 2 and Appendix Figure 6), with five more regressions passed for NDVI, and three more regressions passed for NDWI, clear trends are revealed for both indices and it can be determined that FCH4 tended to correlate positively with NDVI, and negatively with NDWI. 41 Figure 15 US-Myb simple linear regression. Scatterplots of simple linear regression analyses for US-Myb over the 2019. Each row is an individual observation, labelled by the Landsat 8 image acquisition date. Highlighted in red are Pearson correlation values of runs whose p-values are larger than 0.1. 42 Table 4 US-Myb r-value data table. Statistically significant Pearson correlation coefficient values for US-Myb assessed at the α = 0.1 significance level. Positive r-values are highlighted blue, and negative r-values are highlighted orange. Observations where all correlations were not significant are omitted. US-Myb Observation NDVI NDWI NDMI LST 2021-05-08 -0.899 0.837 - -0.873 2021-05-24 - - - -0.852 2021-06-09 - 0.831 - - 2021-06-25 - - 0.952 -0.97 2021-07-11 - - - -0.944 2021-08-28 - - - -0.968 2021-09-13 - -0.853 0.865 - 2021-09-29 0.984 -0.987 - - Figure 16 Boxplots of statistically significant Pearson correlation coefficients for US-Myb. r-value data points are plotted as blue circles, and medians are plotted as the red lines. 43 3.1.4 Simple linear regression – US-WPT Simple linear regression analyses for 2013 at US-WPT, plotted in Figure 17, show that FCH4 was particularly strongly correlated with NDVI and NDWI on June 20th with Pearson r-values of -0.996 and 0.989, respectively. The mean absolute r-values for NDVI, NDWI, NDMI, and LST for US-WPT were 0.927, 0.937, 0.895, and 0.923 respectively. Every observation performed at this site yielded at least one linear regression with a Landsat 8 index that was statistically significant, but more than half of the total correlations were not significant at the α = 0.1 significance level. Summarized in Table 5 and Figure 18, the remaining correlations for NDVI and NDWI were equally split along the r = 0 line, making it difficult to reliably evaluate their relationships with FCH4. Unlike US-Myb, none of the NDVI and NDWI regressions were in the marginal p = 0.1 – 0.15 range. NDMI overall tended towards a negative relationship with FCH4, and LST was consistently positively correlated with FCH4. 44 Figure 17 US-WPT simple linear regression. Scatterplots of simple linear regression analyses for US-WPT over the 2013. Each row is an individual observation, labelled by the Landsat 8 image acquisition date. Highlighted in red are Pearson correlation values of runs whose p-values are larger than 0.1. 45 Table 5 US-WPT r-value data table. Statistically significant Pearson correlation coefficient values for US-WPT assessed at the α = 0.1 significance level. Positive r-values are highlighted blue, and negative r-values are highlighted orange. US-WPT Observation NDVI NDWI NDMI LST 2013-05-26 0.858 -0.887 - - 2013-06-20 -0.996 0.989 - 0.931 2013-07-13 - - -0.871 0.915 2013-08-23 - - 0.864 - 2013-09-24 - - -0.949 - Figure 18 Boxplots of statistically significant Pearson correlation coefficients for US-WPT. r-value data points are plotted as blue circles, and medians are plotted as the red lines. 3.1.5 Summary of simple linear regression results The median Pearson r-values for each Landsat 8 index at all sites are compiled in Figure 19. 46 Figure 19 Median r-value heat map. A summary of median Pearson r-value results from the linear regression analyses plotted as a heat map. For MBPPW1 and MBPPW2, which had multi-year observations from 2021 and 2022, additional median values were calculated which were inclusive of both years. Broadly looking at the aggregated median values for MBPPW1, NDVI and NDMI were negatively related with FCH4, NDWI was positively related with FCH4, and LST’s correlation was more variable between years. For MBPPW2, NDVI and LST were negatively correlated with FCH4, while NDWI and NDMI were positively correlated with FCH4. Landsat 8 index relationships at MBPPW2 were notably more consistent than the other sites. US-Myb had inconclusive relationships for NDVI and NDWI, positive correlation for NDMI, and a negative correlation for 47 LST. However, at the α = 0.15 significance level, US-Myb shows positive correlation trends for NDVI, and negative correlation trends for NDWI. US-WPT also had inconclusive relationships FCH4 for NDVI and NDWI, but a positive correlation for NDMI, and a negative correlation for LST. Further interpretations of these results to characterize the CH4 dynamics at each site are made in the Discussion section. 3.2 FCH4 Prediction with Multiple Linear Regression The predictor variables for the two MLR models are as listed below: Independent variables Simple Model - NDVI - NDWI - NDMI - LST - month.f Stepwise Model - NDVI - NDWI - NDMI - month.f - month.f : NDVI - month.f : NDWI The Stepwise model, built with the criterion of improving the adjusted R2, consisted of the same variables as the Simple model, but with LST excluded, and with the inclusion of interactions between month.f and NDVI and month.f and NDWI. Shapiro-Wilk Tests yielded p-values of 0.884 for the Simple model, and 0.083 for the Stepwise model. With both p-values larger than the 0.05, both model were determined to meet the ordinary least squares normality assumption with residuals following a normal distribution. Predicted-y and residual plots were used to assess the linearity and residual homoscedasticity assumptions of MLR (Figure 20). 48 Figure 20 Predicted-y and residual plots. The top row of plots shows the predicted-y plot (A) and residual plot (B) for the Simple model, and the bottom row shows the predicted-y plot (C) and residual plot (D) of the Stepwise model. The red line in the predicted-y plots represents the 1:1 line, and the red line in the Residual plots represents the zero-line. Data are modelled after MBPPW2 data. The Simple model’s predicted FCH4 values were able to follow the linear 1:1 line with a coefficient of determination (R2) of 0.754, and the Stepwise model’s values were able to follow the 1:1 line with a slightly higher R2 of 0.806 (Table 6). The linearity assumption was concluded to be met based on visual assessment of these plots. The residual plots for both Models showed minor deviations in residual variances from homoscedasticity. Residual variance on the positive half (upper halves of Figure 20 B and D) had a slight fan-pattern in the low range of predicted FCH4 from 0 – 20 nmol/m2/s, but became consistent across the rest of the range from 20 – 100 nmol/m2/s. On the negative half of residual variances (lower halves of Figure 20 B and D), the increasing-spread, fan pattern was more pronounced across most of the y range. Tests were done to incorporate a complex error structure with a weighted variance function to the models, but this caused the models to notably violate the MLR linearity and normality assumptions. Since the 49 residual variances only deviated slightly from equivariance, the models were left un-corrected for heteroscedasticity. Comparing the two models using adjusted R2 to account for their different number of variable terms, the Stepwise model had a higher adjusted R2 of 0.731, while the Simple model had an adjusted R2 of 0.711 (Table 6), indicating that the Stepwise model, with its higher complexity taken into account, produced a better goodness of fit than the Simple model. Both models were trained on the MBPPW2 clustered data from the following days: 2021: 07/14; 08/08; 08/31; 09/25, and 2022: 05/23; 06/15; 07/01; 07/17; 07/26; 09/03; 09/28 (model coefficients listed in Table 6). Table 6 MLR model coefficients and statistics Simple model Stepwise model β0 βNDVI βNDWI βNDMI βLST β0 βNDVI βNDWI ΒNDMI May 44.84 21.81 66.72 167.97 -0.99 11.76 15.13 17.36 124.82 June -19.94 21.81 66.72 167.97 -0.99 1.32 15.92 61.17 124.82 July -15.31 21.81 66.72 167.97 -0.99 26.50 306.93 420.85 124.82 Aug 7.14 21.81 66.72 167.97 -0.99 -173.21 1830.68 1697.78 124.82 Sept -11.87 21.81 66.72 167.97 -0.99 24.89 -32.76 39.90 124.82 R2 0.754 0.806 Adj R2 0.711 0.731 50 MBPPW2 clustered data from 2021-09-09 and 2022-08-18 were used as the test dataset. Evaluating the models on the test dataset, the Simple model was able to predict with an R2 of 0.689 and 0.304, respectively, for the two dates, and the Stepwise model predicted FCH4 with an R2 of 0.741 and 0.960, respectively (Table 7). Both models were well-fit with the 2021-09-09 data. For 2022-08-19, model performance was lower for the Simple model, as indicated by its R2 value 0.304, while the Stepwise model was able to fit the 2022-08-19 data with an R2 of 0.960. For the regression fits, the regression slopes for both models were similar, and less than 1 for 2021-09-09 (Figure 21A), leading them to both over-estimate FCH4 in the mid-to-low range of values, and under-estimate FCH4 in the highest range (> 80 nmol/m2/s) . The Simple model appeared to predict the data slightly better with a regression slope of 0.514, which is slightly higher than the Stepwise model’s slope of 0.469. Figure 21 Predicted FCH4 plots. Predicted FCH4 for 2021-09-09 are plotted in the left plot (A), and predictions for 2022-08-18 are plotted in the right plot (B). The blue and black lines are linear regression lines for the Stepwise model and Simple model respectively, and the red line is the 1:1 line. 51 For 2022-08-18 (Figure 21B), the Simple model again had a regression slope much lower than one, as well as a large positive bias, leading to an over-estimation of FCH4 across the range of values (Table 7). The Stepwise model produced a much better slope approximating 1, but an extremely large positive bias also led it to over-estimated values across the range of FCH4. Based on the 2021-09-09 results, both models performed similarly, with the Stepwise model producing marginally better results with slightly lower bias. Based on the 2022-08-19 results, the Simple model would perform better in the application of ecosystem upscaling due to its lower error bias, and the Stepwise model with its better regression fit and accuracy, would perform better in characterizing and predicting FCH4 spatial distribution. Table 7 MLR model prediction results. R2 Regression slope Regression offset Simple model 2021-09-09 0.69 0.51 43.00 2022-08-18 0.30 0.26 68.66 Stepwise model 2021-09-09 0.74 0.47 39.46 2022-08-18 0.96 1.34 73.32 52 Chapter 4: Discussion The objectives of this study were to investigate whether relationships between remote sensing indices of a wetland surface and the spatial variability of its CH4 emissions could be observed, and then to explore the potential of extrapolating these relationships to estimate FCH4 beyond an EC flux tower’s footprint. The pairing of remote sensing with Footprint-Weighted Flux Maps to achieve these objectives is a novel approach to studying wetland FCH4 spatial dynamics that enables the characterization of FCH4 spatial variability within a single EC flux tower’s footprint. Results from the simple linear regression analyses were used to investigate the first objective by using Landsat 8 indices to characterize the different land surface conditions as they relate to the different FCH4 features identified by the Footprint-weighted Flux Maps. MLR modelling was then used to test the ability to use these indices to predict FCH4 from remote sensing data alone at MBPPW2. 4.1.1 Methane fluxes at MBPPW1 The linear regression results are presented graphically as a heat map in Figure 19. The findings from linear regression analyses at MBPPW1 reveal that this study’s methods may have limited application at wetland sites with particularly low FCH4 as discussed below. In 2021, FCH4 was found to positively correlate with NDVI and LST, while NDWI was negatively correlated. NDMI had an equal split of observations between negative and positive correlations. The signs of the NDVI, NDWI, and LST relationships flipped in 2022, and NDVI, LST, became negatively correlated, while NDWI became positively correlated. Generally speaking, from the full 2021-2022 dataset, higher FCH4 were found from the high NDWI open 53 water areas, but the inconsistencies in the relationships across observations make it difficult to come to a reliable conclusion about the spatial indicators of CH4 emissions at MBPPW1. Explanations for the high inconsistency could be owed to low overall CH4 activity at the site (Figure 5a), or potential differences in CH4 activity across years. MBPPW1 is characterized by low CH4 emissions resulting from high water sulfate levels (Soued et al., in preparation). With sulfate suppressing CH4 production, CH4 hotspots are likely rare occurrences at this site, and without distinct FCH4 features like hotspots to tie CH4 emissions to specific surface conditions, detection of CH4 activity would likely become spatially inconsistent across observations. The major controls of CH4 activity could also evolve inter-annually. With the complexity and temporal variability of CH4 drivers, it’s possible that higher CH4 emissions indeed existed over greener emergent vegetation in 2021, but then shifts in biotic and abiotic drivers in 2022 could have led to higher emissions over the more inundated open water areas of the flux footprint instead (Delwiche et al., 2021; Knox et al., 2021). Box plots of FCH4 over the growing season from both years show moderate annual variation (Figure 22). A t test on inter-annual FCH4 supports this with a p-value of 0.0015, indicating that mean growing season FCH4 for the two years were statistically different. 54 Figure 22 Box plot of FCH4 at MBPPW1. FCH4 data are split into months and years, with pink bars representing 2021 data, and blue bars representing 2022 data. Orange lines within the bars represent the median values. Biotic drivers in particular may be a strong control of inter-annual FCH4 variability at this site. MBPPW1 is dominated by Typha around the flux tower site (Figure 1), and from the random forest model output of CH4 driver importance (Figure 6), the strongest predictors of FCH4 were H, NEE, and TS 10 cm beneath the sediment-water interface. Typha plant biomass over the growing season is highly correlated with NEE in temperate marshes, and H and TS can be strongly controlled by canopy density (Bonneville et al., 2008; Tanaka & Hashimoto, 2006). Any spatio-temporal variability in Typha population could impact their contribution to ecosystem productivity and canopy biomass, and ultimately influence the spatio-temporal variability of CH4, as detected in this linear regression analysis. 55 MBPPW1 was established in 2021, so only two years of growing season data were available at the time of writing this thesis. Further analysis at this site over additional years could help explain more of the inter-annual variability of surface conditions and FCH4. 4.1.2 Methane fluxes at MBPPW2 Median Pearson r-values representing Landsat-FCH4 relationships for NDVI, NDWI, and NDMI at MBPPW2 agreed consistently across the two years of study for MBPPW2. FCH4 remained negatively correlated with NDVI, and positively correlated with NDWI and NDMI (Figure 19). These results indicate that at this site, higher CH4 emissions were typically found over the open water where NDVI is low and NDWI is high. The positive relationship found with NDMI indicate that within the vegetated areas, higher CH4 emissions were found with higher soil saturation and vegetation moisture. LST was inconsistent, which was likely due to seasonal variations in emergent vegetation and aquatic vegetation within the open water. Further discussions on the reliability of LST are made in the linear regression limitations section. CH4 emissions being higher over open water rather than from the surrounding emergent vegetation was an unexpected result considering the well-understood role of plants in the promotion of CH4 production through carbon substrate production, and facilitation of CH4 release as an emission pathway though aerenchymous plant tissue (Colmer, 2003; Joabsson et al., 1999; Knox et al., 2021; McNicol et al., 2017; Schipper & Reddy, 1994). Many studies in wetlands with heterogeneous mixes of open water and emergent vegetation often found CH4 emissions from the open water to be lower relative to emissions from the emergent vegetation (Li et al., 2020; McNicol et al., 2017; Strachan et al., 2015). Strachan et al. (2015) found open 56 water CH4 fluxes in a temperate freshwater marsh in the range of 0.5 ± 0.4 g C m-2 d-1 as compared to CH4 fluxes from emergent Typha in the range of 1.0 ± 1.2 g C m-2 d-1. Li et al. (2020) calculated open water CH4 fluxes from an urban wetland to be 1.73 g C-CO2 m-2 yr-1, and vegetated fluxes to be 4.26 73 g C-CO2 m-2 yr-1. And at US-Myb, McNicol et al. (2017) found open water fluxes in the range of 2.9 ± 0.5 g C m-2 yr-1, and vegetated fluxes in the range of 94.8 ± 12. g C m-2 yr-1. At MBPPW2, different wetland processes could be occurring to either inhibit vegetated CH4 emissions or promote open water emissions. The emergent vegetation themselves could be inhibiting CH4 emissions. While the aerenchymous plant tissue in Typha provide vertical transport pathways for CH4, their roots can also simultaneously introduce oxygen into the surrounding anaerobic soil layer during respiration (Watson et al., 1997). The supplied oxygen can stimulate methanotrophs around the roots to increases CH4 oxidation rates, which offsets CH4 production rates, and reduces overall emissions from the emergent vegetation (Joabsson et al., 1999; Watson et al., 1997). Additional research in the future can be performed using stable isotope analysis to investigate CH4 oxidation rates between open water and emergent vegetation. Open water CH4 emissions could also be elevated due to the contributions from sub-surface macrophytes. During the warmer months, large mats of macrophytes will develop near the open water surface at MBPPW2. These mats can be seen from high resolution drone footage (Figure 23) in the open water to the south of the flux tower. This aquatic vegetation could be supplying carbon substrates to stimulate CH4 production in the open water soils, as well as providing plant mediated vertical transport of CH4 through the water column. If the rates of these processes over 57 aquatic conditions are comparable to those over the emergent vegetation, the added contribution from WTD as a driver of FCH4, and the varying periods and extents of soil inundation between the open water and emergent vegetation areas could explain the higher emissions from the open water. Since these macrophyte mats are located beneath the water surface, the NDVI calculated from these open water pixels would be lower than from the surrounding emergent vegetation, leading to the negative NDVI relationships found with FCH4. Figure 23 Drone photography of MBPPW2. Drone footage was taken on 2022/09/29 by Paige Kowal from the Ducks Unlimited Institute of Wetland and Waterfowl Research. The red star marks the location of the flux tower. MBPPW2’s linear regression results were also supported by an additional flux footprint analysis done with LI-COR Biosciences’ Tovi Flux Footprint Allocation (FFA) tool. FFA allows users to graphically designate regions of interest (ROI), and allocate flux data based on the percentage overlap between flux footprints and ROI area. Performed on MBPPW2’s 2021 EC data, the 58 wetland’s surface was separated into open water and emergent vegetation ROI (bottom left of Figure 24). The graph of average diurnal ROI-allocated-FCH4 in the bottom-middle of Figure 24 showed that across a 24-hour cycle, the average CH4 emissions from the open water were higher than the average emission from emergent vegetation. This supports the findings of higher FCH4 from the open water from this study. To corroborate these findings, future research could further explore the spatial variability in FCH4 by conducting chamber measurements over the emergent vegetation and open water areas. Figure 24 Tovi Flux Footprint Allocation (FFA) results for FCH4 at MBPPW2. FFA analysis was performed using 2021 data from MBPPW2. EC data cleaning and gap-filling were performed within the Tovi workflow. The land surface is separated into open water in blue, and emergent vegetation in yellow. The top plot is a CH4 time series with data points color-coded to the ROI that they are allocated into. Bottom left image shows the ROI allocation areas. The bottom middle plot shows the average diurnal CH4 fluxes for both ROI-allocated datasets. Since we don’t yet have flux chamber data available, a parallel FFA analysis was conducted with CO2 fluxes (Figure 25) to help support the FCH4-FFA findings from Figure 24. This run of the analysis tool successfully identified stronger drawdown of CO2 over the emergent vegetation ROI during daytime from 12pm and 7pm, as should be expected, highlighting the efficacy of 59 FFA as a tool. It additionally observed net CO2 uptake over the open water ROI during the daytime, suggesting the presence of aquatic vegetation in MBPPW2’s open water. This aquatic vegetation could be the same macrophytes seen in Figure 23, and their contribution to GPP processes lends confidence to the plausibility of their role in facilitating strong open water FCH4. Figure 25 Tovi FFA results for CO2 at MBPPW2. Footprint-allocated data plotted are CO2 fluxes. 4.1.3 Methane fluxes at US-Myb Although there were some discrepancies in the linear regression analysis results at the 10% significance level, the overall results were in line with previous CH4 research conducted at US-Myb (McNicol et al., 2017). Across observations, NDVI and NDWI r-values were equally split into negative and positive relationship ends (Figure 16). NDMI and LST, however, were more consistent, with median Pearson r-values being positively correlated for NDMI, and negatively correlated for LST (Figure 19). These results align well with this site’s literature describing the site’s FCH4 spatial structures. Flux chamber studies by McNicol et al. (2017) had discovered far 60 stronger CH4 emissions coming from typha-populated emergent vegetation than from the open water (94.8 ± 12.3 vs. 2.9 ± 0.5 g C m-2 yr-1). Within the emergent vegetation, the positive NDMI correlation relates high FCH4 with healthy plants as indicated by plant tissue moisture, which supports the well-understood processes of plant productivity facilitating CH4 emissions through plant-mediated transport and methanogen activity through supply of carbon substrate into the soil (Gao, 1996; Joabsson et al., 1999; Oertel et al., 2016; Rokni et al., 2014). Regarding NDVI and NDWI, while their results were inconclusive at the 10% significance level, loosening the statistical significance to α = 0.15 included 8 more significant correlations to paint a much more consistent image of Landsat-FCH4 relationships (Appendix Figure 6). From these values, NDVI was clearly related positively to FCH4, and NDWI was negatively correlated. Again, these results also align well with the findings from McNicol et al. (2017), with both agreeing on stronger CH4 emissions coming from greener, emergent vegetation as compared to those coming from the open water. 15% is not a commonly used significance level, but assessing the results at α = 0.15 was able to provide additional support to this, and previous studies’ findings of higher FCH4 over emergent vegetation. As a note on the agreement between this and previous study’s results, McNicol et al. (2017) based their study on data from 2014, while this study used data from 2021. The time elapsed between studies is not expected to lead to major differences in the site’s CH4 dynamics, but it is worth noting that different periods were being investigated in these two studies. The reason for selecting 2021 in this thesis was that 2021 provided the largest number of statistically significant 61 observations. For 2014, only 3 out of 40 regressions would have been significant at the α = 0.1 significance level (Appendix Figure 7). A possible explanation for why US-Myb observations in particular struggled with significance is that this site exhibited strong wind direction bias. Plotting the wind direction data from this site, and using MBPPW2 as a basis of comparison (Figure 26), a large majority of the EC flux measurements were measured when the wind was from a westerly direction at around 270˚ wind direction. Figure 26 Wind direction plot. Half-hourly EC wind direction data plotted for US-Myb on the left, and MBPPW2 on the right for comparison from their respective study period(s). Other sites in this study, like MBPPW2, had a much more even sampling of wind directions. This strong bias at US-Myb caused Footprint-weighted Flux Maps to model narrow footprints with poor spatial coverage of the land surface (example Flux Map in Figure 27), that may not be conducive to the 30-m spatial resolution coarsening and 5 cluster number grouping done in this study. This also points to a limitation of the clustered-linear regression method and suggests that upscaling may be more reliable at sites with winds coming from a broader range of directions. 62 Figure 27 Footprint-weighted Flux Map output. Flux Map for US-Myb was generated for the period from 2021-05-26 to 2021-06-23. 4.1.4 Methane fluxes at US-WPT 2013 was the only year with available data from US-WPT that overlapped with the Landsat 8 satellite operational dates. Similar to US-Myb, linear regression results for NDVI and NDWI were split equally in positive and negative r-value signs, so their spatial pattern trends could not reliably determined (Figure 18). However, NDMI was consistent in its negative correlation and LST was consistent in its positive correlation. At US-WPT, NDMI had a particularly strong (negative) correlation with NDWI (Figure 28). Figure 28 NDMI vs, NDWI scatter plot. Plotted data are clustered Landsat 8 NDMI and NDWI values from US-WPT. 63 The negative NDMI relationship at this site could indicate that, similar to MBPPW2, stronger CH4 emissions were coming from the open water rather than from emergent vegetation. According to previous studies performed at this site, US-WPT had large aquatic macrophyte populations in its open water, and exceptionally strong CH4 emissions coming from the aquatic water-lilies and macrophytes (Chu et al., 2014, 2016). Chu et al. (2014)’s study period at US-WPT spanned from 2011 to 2013, so their findings are relevant, and help support the interpretation of stronger open water FCH4 at this site. 4.1.5 Linear regression analysis – key findings The spatially clustered linear regression analyses in this study were able to demonstrate repeatable relationships coupling remotely sensed Landsat 8 images with FCH4 spatial dynamics across multiple sites. Agreement between these findings and results from either previous studies or supplemental analyses supports the validity of this novel approach and the potential for its application in characterizing within-footprint FCH4 spatial structures, and beyond that, ecosystem upscaling and FCH4 prediction. While this approach has proven effective in investigating the first objective in this study, further improvements must be made to address its limitations and expand it to a broader range of sites in order to improve its reliability and wider application. 4.1.5a Linear regression analysis - statistical limitations and future work The use of hierarchical clustering to reduce the observation sample sizes down to 5 data points per regression makes it challenging to assess the statistical assumptions of simple linear regression. One strategy that could have helped address this limitation was to use Spearman rank 64 correlation instead of Pearson correlation. Spearman rank correlation is a nonparametric measure that is not limited by the assumptions of normality and homoscedasticity, and is typically a better metric than Pearson correlation for small sample sizes (Kendall & Stuart, 1973). However, because Spearman is a rank correlation, its correlation coefficients do not represent the linear relationships between variables, but rather the monotonic ranks of the data. For the purposes of this study, Pearson presented a more useful measure of relationship type and strength between Landsat 8 vs. FCH4. Spearman correlation values were still checked in each regression analysis to make sure there were no major departures from the tested Pearson values. The choice of a 0.10 significance level over a more commonly used 0.05 significance level may also have impacted the reliability of these results. The reason for choosing the looser standard was simply to increase the number of sample points made passable to the final analyses. This can be compensated for in future adjustments to this study by increasing the temporal sampling frequency. A possible way to achieve this is by incorporating additional remote sensing data from the newly launched Landsat 9 satellite. Landsat 9 has the same instruments, bands, and remote sensing products as Landsat 8, but orbits in an offset position, allowing it to capture images of a site on different days from Landsat 8. Incorporating both satellite products will produce high temporal sampling frequencies, larger datasets, and possibly enable the use of the 0.05 significance level to create more reliable results. 4.1.5b Linear regression analysis – practical limitations and future work Many of the NDWI calculations came back as negative values over the inundated and open water areas of the wetland sites. This can arise from presence of organic matter in the water, or mixture 65 of water with wet soils (Singh et al., 2015). While the relative values of the pixels were more of the focus in this study, inaccuracies from improperly sensed water features could cascade across calculations and influence the results. This study could be improved in the future by incorporating a Modified Normalized Difference Water Index, calculated from the green and shortwave infrared bands rather than the green and near-infrared bands, which has been argued to be more suitable for water features that are closely surrounded by land formations (Singh et al., 2015; Xu, 2006). LST results were highly variable across observations for many of the sites. This can be owed to fact that in heavily vegetated sites with dense canopy layers, surface temperature detected from the top of the canopy may not be representative of TS beneath the canopy floor. Further, variability in macrophyte abundance in the open waters may influence broader temperature dynamics between emergent vegetation and open water. Periods with densely-populated, near-surface vegetation in the open water may lead to higher LST over the open water than emergent vegetation, but conversely, when there is sparse aquatic vegetation, the open water, with its higher heat capacity and thermal conductivity, may have a lower LST than the emergent vegetation. This variability could have been what influenced the linear regression results, to make LST-FCH4 relationships inconsistent across observations. LST could still be a strong tool in this research, but it may not be the most dominant predictor for within-site spatial variability in FCH4. The Footprint-weighted Flux Maps have the capability to analyze footprint-weighted sensible and latent heat flux, and could be a good starting point to investigate spatial variability of LST in wetlands. 66 The final limitation is the wide range of time-scales in which FCH4 interacts with its drivers. FCH4 can interact with drivers on a seasonal scale where it lags TS by 5 ± 15 days, lags water table dynamics by 17 ± 11 days, and lags GPP by 12.5 ± 23 days (Knox et al., 2021). On a multi-day scale, FCH4 is near-synchronous with TS, lags WTD by ~8 hours, and lags GPP by 1 – 4 hours (Knox et al., 2021). The variability within and across each variable makes it difficult to accurately match the timing between FCH4 measurements and Landsat 8 detection of surface conditions. This study tried to be as inclusive of all the interaction windows as possible with its 4-week Flux Map window, but this may be too rigid of a structure. More accurate observations could be made possible if each FCH4 driver, with its represented Landsat 8 index, were matched with a separate Footprint-weighted Flux Map created from a fitted data-window that captures the appropriate time-scale of the driver’s interaction with FCH4. The current state of this study’s clustered linear regression analysis lays a strong foundation for its novel applications in spatial FCH4 characterization, spatial upscaling and FCH4 prediction. Continuing to build on this foundation by incorporating the future work mentioned above can improve the study’s reliability, allowing it to become a useful tool for broader applications to future research in freshwater wetland FCH4. 4.2 MLR modelling The MLR modelling in this study addresses the second objective of this thesis by exploring the application of spatially clustered Landsat-FCH4 relationships in ecosystem upscaling and prediction. The two MLR models tested in this study were able capture a considerable amount of variation in the data, and their results present a promising framework for the development of a 67 novel approach for upscaling FCH4 dynamics in wetland ecosystems using only remote sensing and a single EC flux tower. The reliability of these models, however, still require improvements to address a number of limitations, but before discussing these limitations, an alternative approach to CH4 prediction was first tested to examine a functional next step in the application of clustered MLR modelling. 4.2.1 Applicability of MLR predictions MLR on the test dataset produced fairly accurate results, but this study’s iteration of the method represented just the first step in predicting CH4 beyond the range of EC flux towers. Currently, the model is being built and tested on data that were hierarchical clustered based on information from the Footprint-weighted Flux Maps, which are inextricably linked with EC flux towers. Applying these current models to prediction using raw, un-clustered remote sensing data yielded poor model performances (Figure 29). Pearson r-values were in the range of 0.254 to 0.448, and predicted FCH4 clearly did not follow a linear 1:1 relationship with observed FCH4. 68 Figure 29 MLR prediction on un-clustered data Plotted FCH4 are un-clustered, but resolution- and footprint-matched on the two test dates at MBPPW2. Black dotted lines are 1:1 lines, and solid orange lines are linear regression fit lines. These results could be improved with an alteration on this study’s clustering method by reversing the roles of the two datasets to instead perform hierarchical clustering on Landsat 8 data information, and then having the cluster labels be applied to Footprint-weighted Flux Maps. Preliminary attempts at this method were briefly explored by simply using the existing MLR models built in this study, and predicting on Landsat 8 indices clustered on NDVI data. The resulting clusters here would then be representative of NDVI features in the MBPPW2 footprint coverage (Figure 30). 69 Figure 30 MBPPW2 Landsat 8 NDVI images Image for 2021-09-09 is on the left, and 2022-08-18 is on the right. The red circle marks the approximate flux footprint climatology area. The prediction results on the NDVI-clustered dataset (Figure 31) were, expectedly, far less accurate than the results from the original Flux-Map-clustered dataset (Figure 21), but improved on the raw, un-clustered dataset (Figure 29). There was an issue where the NDVI measurements within the MBPPW2 flux footprint on 2021-09-09 (Figure 30) were much more strongly contrasted than NDVI on 2022-08-18, and this appeared to have affected the efficacy of clustering on 2022-08-18. The prediction results on the clustered data, informed by the much more homogenous 2022-08-18 NDVI, produced poor predictions for 2022-08-08 (Figure 31B), with R2 of 0.02 and 0.009 for the Simple and Stepwise models respectively. The 2021-09-09 NDVI-informed predictions (Figure 31A) appeared much better with improved goodness of fits (R2 of 0.926 for Main effects model and 0.930 for Reduced model), but both models produced regression slopes that were too steep, leading to vastly overestimated FCH4 on the higher end. 70 Figure 31 MLR prediction results on NDVI-clustered data. Predictions for 2021-09-09 are plotted on the left (A) and predictions for 2022-08-18 are plotted on the right (B). The red line represents the 1:1 line. Simple model results are plotted in black and Stepwise model results are plotted in blue. While further research is needed, this approach is worth exploring further as it is not reliant on EC data to pre-process any of the predictor data. This would freely allow the models to predict in wetland areas without requiring EC data. Future improvements to this approach would entail deeper investigations into more effective ways to inform hierarchical clustering with Landsat 8 data, such as using different Landsat indices under different conditions, and to build new MLR models from the ground up using the new clustering approach. More importantly, however, overall improvements may be made by first addressing the limitations in the original MLR models. 4.2.2 Limitations of MLR - temporal dependency To obtain enough datapoints to populate the training datasets, data had to be aggregated over multiple observations in time. However, FCH4 and many of the predictor variables display strong seasonal and annual temporal patterns. Including month-of-the-year as a predictor helped 71 mitigate the influence of temporal dynamics, but additional steps could be explored to make further improvements. The training dataset sampled repeated measurements from the same spatial surface within the flux footprint over irregular time intervals, potentially introducing temporal autocorrelation to the dataset. ACF tests confirmed this as a problem in both MLR models. To account for this, a generalized least squares model could be used to incorporate an AR(1) correlation structure in the models. Specifically, in the R nlme package’s gls function, the continuous AR(1) correlation structure, corCAR1, could be added to properly account for irregular sampling arising from limitations in remote sensing (Pinheiro et al., 2019). Additionally, temporal de-trending of Landsat 8 and Flux Map data between observations with the use of a high-pass Fourier Transform could potentially filter out much of the temporal signal from the dataset. The inclusion of Year as an additional variable could also help accommodate for inter-annual variation in the data. It was attempted to include Year in the models, but was not successfully implemented due to insufficient coverage of all growing season months within each year. Gathering enough data to allow Year to be included could help reduce the stronger error biases in the 2022-08-18 predictions. 4.2.3 Limitations of MLR - multicollinearity and spatial dependency There was strong multicollinearity between the predictor variables. From a visual assessment of a scatterplot matrix (Figure 32), NDVI was particularly strongly correlated with both NDWI and NDMI, while NDWI and NDMI were moderately correlated. This was supported by Variance 72 Inflation Factor (VIF) tests on the two tested models, where VIF values less than 1 indicated no collinearity, between 1 – 5 indicated moderate collinearity, and greater than 5 indicated high collinearity. For the Main Effects model, month.f, NDMI, and temp had VIF values between 1 and 5 (moderately collinear), and NDVI and NDWI had VIF values above 5 (strongly collinear). The Reduced model exhibited much stronger collinearity, with each of its variables except NDMI having a VIF value in the range of 9.2 – 167.4. This study used all relevant indices to investigate FCH4 and its drivers and thus did not address multicollinearity in its models. Future improvements could be made to these methods to adjust for multicollinearity by exploring Principal Component Analysis to transform the predictor variables, or building alternate MLR models that exclude highly collinear variables altogether. 73 Figure 32 MBPPW2 Scatterplot matrix of MLR variables Plotted data are from the clustered, and statistically significant 2021 and 2022 datasets MBPPW2. Spatial dependency is another factor that could have influenced the reliability of the MLR models, as the sampled data all came from neighbouring pixels within single flux footprint areas. Including the corSpatial correlation structure in the MLR models could have addressed spatial autocorrelation, but it was not done so in favour of maintaining low model complexity. Instead, the hierarchical clustering approach was relied upon to ensure that neighbouring clusters of data were sufficiently distinct from each other to mitigate the spatial correlation between data points. 74 4.2.4 Additional modelling options This study used all-fixed-effect MLR models, but alternatively, linear mixed-effect MLR models with month.f and its interactions as random effects may produce better results. However, committing to only using rigid MLR models may limit the robustness of this study. As noted in the MLR Methods section, the heteroscedasticity assumption was moderately violated, and the typical work-around by including a weighted variance function created even more severe linearity problems (Figure 33). Figure 33 Predicted-y plot with weighted variance model Sticking with MLR models may therefore not be an effective approach. Other options such as random forest models or multilayer perceptron artificial neural networks do not have such assumption demands, and may provide much-needed model flexibility for this data. 75 Chapter 5: Conclusion Linear regression analysis on spatially clustered Landsat 8 and Footprint-weighted Flux Map data was able to find strong correlations relating remotely sensed surface conditions (NDVI, NDWI, NDMI, and LST) to FCH4. While observations were not all statistically significant at the 10% level, taking the median Pearson r-values across observations was able to produce meaningful results. This study was able to characterize stronger open water CH4 emissions at the MBPPW2 and US-WPT wetland sites, and conversely, stronger emergent vegetation CH4 emissions at the US-Myb site. It was possible to weakly detect stronger CH4 emissions over open water at MBPPW1. These characterizations were able to be supported by either existing literature, as was the case for US-Myb and US-WPT, or supplemental analyses, as was the case for MBPPW2. The findings from all four sites also demonstrate the strengths and limitations of the broad application of this study’s methods. For sites like MBPPW2 and US-WPT with distinct CH4 emission sources and adequate footprint coverage of the land surface, strong characterizations of spatial distributions of FCH4 can be obtained. For sites with low CH4 emissions and strongly biased wind directions, like MBPPW1 and US-Myb, respectively, it is difficult to produce consistent, and statistically significant clustered-regression correlations. In future applications of these methods for any of FCH4 characterization, upscaling, or prediction, their strengths and limitations should taken into account to assess the suitability of a wetland site before proceeding. The FCH4 modelling done in this study is the first of its type to utilize within-footprint FCH4 variability to enable spatial estimates of FCH4 at relatively fine scales. Both the Simple and the 76 Stepwise MLR models were able to predict FCH4 on the two testing dates, 2021-09-09 and 2022-08-18, with high R2, meaning much of the FCH4 variance was captured with either model. However, despite the high R2 value, both models had issues with accuracy. Both models exhibited positive bias for both test dates, and both of their regression slopes were too small for 2021-09-09, with the Simple model continuing to produce a regression slope < 1 for 2022-08-08. The Stepwise model predicted 2022-08-08 FCH4 much more precisely with a regression slope approximating 1, but had a much stronger bias. Overall, the Stepwise model demonstrated that it would perform better with the objective of characterizing FCH4 spatial structures due to its better-fit-of-data enabling better spatial prediction of relative fluxes. The Simple model, on the other hand, with its lower positive bias and lesser likelihood to overestimate FCH4, might perform better in the application to ecosystem upscaling. 5.1 Research objectives and next steps The results of this research were able to supported the first objective, which was to investigate whether significant relationships between Landsat 8 remote sensing indices and FCH4 were detectable. Although consistency remains an issue, strongly correlated and statistically significant relationships were demonstrated for each site. Regarding the second objective, the study showed a strong potential of extrapolating these relationships to estimating FCH4 beyond the EC flux towers, but considerable work remains to be done. However, the progress made in this study has built a strong foundation for this objective, and there are very promising future steps that can be taken to continue to build on this established framework. This includes incorporating Landsat 9 into the dataset to increase temporal sampling frequency, adapting this study’s Flux Map-informed clustering method to instead cluster based on Landsat data, 77 designing flexible data-windows for the Flux Map to accommodate each predictor’s particular FCH4 relationship timescale, and exploring less-rigid predictive models like random forest and multilayer perceptron neural networks. In addition, further work can be done to improve the reliability of this study’s methods by improving NDWI accuracy with a Modified-NDWI calculation, and by improving the implementation of LST through a deeper understanding of its spatial dynamics in wetlands by applying additional Footprint-weighted Flux Map analyses to H and LE data. With enough fine-tuning of this study’s methods, the scope of their application may be widened to include more wetland sites across various wetland types. Improving the broader application of this study to being able to reliably estimate FCH4 in under-studied, cost-prohibitive, or difficult to reach wetland areas will have tremendous benefits to advancing our understanding of wetland CH4 dynamics, and the use of wetlands as a nature-based climate solution. 78 References Amiro, B. D. (1998). Footprint climatologies for evapotranspiration in a boreal catchment. Agricultural and Forest Meteorology, 90(3), 195–201. https://doi.org/10.1016/S0168-1923(97)00096-8 Badiou, P., McDougal, R., Pennock, D., & Clark, B. (2011). Greenhouse gas emissions and carbon sequestration potential in restored wetlands of the Canadian prairie pothole region. Wetlands Ecology and Management, 19(3), 237–256. https://doi.org/10.1007/s11273-011-9214-6 Baldocchi, D. (2014). Measuring fluxes of trace gases and energy between ecosystems and the atmosphere—The state and future of the eddy covariance method. Global Change Biology, 20(12), 3600–3609. https://doi.org/10.1111/gcb.12649 Barr, A. G., Richardson, A. D., Hollinger, D. Y., Papale, D., Arain, M. A., Black, T. A., Bohrer, G., Dragoni, D., Fischer, M. L., Gu, L., Law, B. E., Margolis, H. A., McCaughey, J. H., Munger, J. W., Oechel, W., & Schaeffer, K. (2013). Use of change-point detection for friction–velocity threshold evaluation in eddy-covariance studies. Agricultural and Forest Meteorology, 171–172, 31–45. https://doi.org/10.1016/j.agrformet.2012.11.023 Bieniada, A., & Strack, M. (2021). Steady and ebullitive methane fluxes from active, restored and unrestored horticultural peatlands. Ecological Engineering, 169, 106324. https://doi.org/10.1016/j.ecoleng.2021.106324 Bloom, A. A., Palmer, P. I., Fraser, A., Reay, D. S., & Frankenberg, C. (2010). Large-Scale Controls of Methanogenesis Inferred from Methane and Gravity Spaceborne Data. Science, 327(5963), 322–325. https://doi.org/10.1126/science.1175176 79 Bonneville, M.-C., Strachan, I. B., Humphreys, E. R., & Roulet, N. T. (2008). Net ecosystem CO2 exchange in a temperate cattail marsh in relation to biophysical properties. Agricultural and Forest Meteorology, 148(1), 69–81. https://doi.org/10.1016/j.agrformet.2007.09.004 Bridgham, S. D., Cadillo-Quiroz, H., Keller, J. K., & Zhuang, Q. (2013). Methane emissions from wetlands: Biogeochemical, microbial, and modeling perspectives from local to global scales. Global Change Biology, 19(5), 1325–1346. https://doi.org/10.1111/gcb.12131 Chen, B., Black, T. A., Coops, N. C., Hilker, T., (Tony) Trofymow, J. A., & Morgenstern, K. (2009). Assessing Tower Flux Footprint Climatology and Scaling Between Remotely Sensed and Eddy Covariance Measurements. Boundary-Layer Meteorology, 130(2), 137–167. https://doi.org/10.1007/s10546-008-9339-1 Chen, B., Coops, N. C., Fu, D., Margolis, H. A., Amiro, B. D., Black, T. A., Arain, M. A., Barr, A. G., Bourque, C. P.-A., Flanagan, L. B., Lafleur, P. M., McCaughey, J. H., & Wofsy, S. C. (2012). 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, 124, 742–755. https://doi.org/10.1016/j.rse.2012.06.007 Chu, H., Chen, J., Gottgens, J. F., Desai, A. R., Ouyang, Z., & Qian, S. S. (2016). Response and biophysical regulation of carbon dioxide fluxes to climate variability and anomaly in contrasting ecosystems in northwestern Ohio, USA. Agricultural and Forest Meteorology, 220, 50–68. https://doi.org/10.1016/j.agrformet.2016.01.008 80 Chu, H., Chen, J., Gottgens, J. F., Ouyang, Z., John, R., Czajkowski, K., & Becker, R. (2014). Net ecosystem methane and carbon dioxide exchanges in a Lake Erie coastal marsh and a nearby cropland: CH4 and CO2 fluxes in a freshwater marsh. Journal of Geophysical Research: Biogeosciences, 119(5), 722–740. https://doi.org/10.1002/2013JG002520 Chu, H., Luo, X., Ouyang, Z., Chan, W. S., Dengel, S., Biraud, S. C., Torn, M. S., Metzger, S., Kumar, J., Arain, M. A., Arkebauer, T. J., Baldocchi, D., Bernacchi, C., Billesbach, D., Black, T. A., Blanken, P. D., Bohrer, G., Bracho, R., Brown, S., … Zona, D. (2021). Representativeness of Eddy-Covariance flux footprints for areas surrounding AmeriFlux sites. Agricultural and Forest Meteorology, 301–302, 108350. https://doi.org/10.1016/j.agrformet.2021.108350 Colmer, T. D. (2003). Long-distance transport of gases in plants: A perspective on internal aeration and radial oxygen loss from roots: Gas transport in plants. Plant, Cell & Environment, 26(1), 17–36. https://doi.org/10.1046/j.1365-3040.2003.00846.x D′Acunha, B., Morillas, L., Black, T. A., Christen, A., & Johnson, M. S. (2019). Net Ecosystem Carbon Balance of a Peat Bog Undergoing Restoration: Integrating CO 2 and CH 4 Fluxes From Eddy Covariance and Aquatic Evasion With DOC Drainage Fluxes. Journal of Geophysical Research: Biogeosciences, 124(4), 884–901. https://doi.org/10.1029/2019JG005123 Delwiche, K. B., Knox, S. H., Malhotra, A., Fluet-Chouinard, E., McNicol, G., Feron, S., Ouyang, Z., Papale, D., Trotta, C., Canfora, E., Cheah, Y.-W., Christianson, D., Alberto, M. C. R., Alekseychik, P., Aurela, M., Baldocchi, D., Bansal, S., Billesbach, D. P., Bohrer, G., … Jackson, R. B. (2021). FLUXNET-CH4: A global, multi-ecosystem dataset 81 and analysis of methane seasonality from freshwater wetlands [Preprint]. Biosphere – Biogeosciences. https://doi.org/10.5194/essd-2020-307 Ermida, S. L., Soares, P., Mantas, V., Göttsche, F.-M., & Trigo, I. F. (2020). Google Earth Engine Open-Source Code for Land Surface Temperature Estimation from the Landsat Series. Remote Sensing, 12(9), 1471. https://doi.org/10.3390/rs12091471 Etminan, M., Myhre, G., Highwood, E. J., & Shine, K. P. (2016). Radiative forcing of carbon dioxide, methane, and nitrous oxide: A significant revision of the methane radiative forcing. Geophysical Research Letters, 43(24). https://doi.org/10.1002/2016GL071930 Forbrich, I., Kutzbach, L., Wille, C., Becker, T., Wu, J., & Wilmking, M. (2011). Cross-evaluation of measurements of peatland methane emissions on microform and ecosystem scales using high-resolution landcover classification and source weight modelling. Agricultural and Forest Meteorology, 151(7), 864–874. https://doi.org/10.1016/j.agrformet.2011.02.006 Gao, B. (1996). NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sensing of Environment, 58(3), 257–266. https://doi.org/10.1016/S0034-4257(96)00067-3 Griscom, B. W., Adams, J., Ellis, P. W., Houghton, R. A., Lomax, G., Miteva, D. A., Schlesinger, W. H., Shoch, D., Siikamäki, J. V., Smith, P., Woodbury, P., Zganjar, C., Blackman, A., Campari, J., Conant, R. T., Delgado, C., Elias, P., Gopalakrishna, T., Hamsik, M. R., … Fargione, J. (2017). Natural climate solutions. Proceedings of the National Academy of Sciences, 114(44), 11645–11650. https://doi.org/10.1073/pnas.1710465114 82 IPCC. (2022). Global Warming of 1.5°C: IPCC Special Report on Impacts of Global Warming of 1.5°C above Pre-industrial Levels in Context of Strengthening Response to Climate Change, Sustainable Development, and Efforts to Eradicate Poverty (1st ed.). Cambridge University Press. https://doi.org/10.1017/9781009157940 Joabsson, A., Christensen, T. R., & Wallén, B. (1999). Vascular plant controls on methane emissions from northern peatforming wetlands. Trends in Ecology & Evolution, 14(10), 385–388. https://doi.org/10.1016/S0169-5347(99)01649-3 Kayranli, B., Scholz, M., Mustafa, A., & Hedmark, Å. (2010). Carbon Storage and Fluxes within Freshwater Wetlands: A Critical Review. Wetlands, 30(1), 111–124. https://doi.org/10.1007/s13157-009-0003-4 Kljun, N., Calanca, P., Rotach, M. W., & Schmid, H. P. (2015). A simple two-dimensional parameterisation for Flux Footprint Prediction (FFP). Geoscientific Model Development, 8(11), 3695–3713. https://doi.org/10.5194/gmd-8-3695-2015 Knox, S. H., Bansal, S., McNicol, G., Schafer, K., Sturtevant, C., Ueyama, M., Valach, A. C., Baldocchi, D., Delwiche, K., Desai, A. R., Euskirchen, E., Liu, J., Lohila, A., Malhotra, A., Melling, L., Riley, W., Runkle, B. R. K., Turner, J., Vargas, R., … Jackson, R. B. (2021). Identifying dominant environmental predictors of freshwater wetland methane fluxes across diurnal to seasonal time scales. Global Change Biology, 27(15), 3582–3604. https://doi.org/10.1111/gcb.15661 Kormann, R., & Meixner, F. X. (2001). An Analytical Footprint Model For Non-Neutral Stratification. Boundary-Layer Meteorology, 99(2), 207–224. https://doi.org/10.1023/A:1018991015119 83 Lee, S., Black, T. A., Nyberg, M., Merkens, M., Nesic, Z., Ng, D., & Knox, S. H. (2021). Biophysical Impacts of Historical Disturbances, Restoration Strategies, and Vegetation Types in a Peatland Ecosystem. Journal of Geophysical Research: Biogeosciences, 126(10). https://doi.org/10.1029/2021JG006532 Lee, S.-C., Christen, A., Black, A. T., Johnson, M. S., Jassal, R. S., Ketler, R., Nesic, Z., & Merkens, M. (2017). Annual greenhouse gas budget for a bog ecosystem undergoing restoration by rewetting. Biogeosciences, 14(11), 2799–2814. https://doi.org/10.5194/bg-14-2799-2017 Li, X., Wahlroos, O., Haapanala, S., Pumpanen, J., Vasander, H., Ojala, A., Vesala, T., & Mammarella, I. (2020). Carbon dioxide and methane fluxes from different surface types in a created urban wetland. Biogeosciences, 17(13), 3409–3425. https://doi.org/10.5194/bg-17-3409-2020 Lu, W., Xiao, J., Liu, F., Zhang, Y., Liu, C., & Lin, G. (2017). Contrasting ecosystem CO 2 fluxes of inland and coastal wetlands: A meta‐analysis of eddy covariance data. Global Change Biology, 23(3), 1180–1198. https://doi.org/10.1111/gcb.13424 McNicol, G., Sturtevant, C. S., Knox, S. H., Dronova, I., Baldocchi, D. D., & Silver, W. L. (2017). Effects of seasonality, transport pathway, and spatial structure on greenhouse gas fluxes in a restored wetland. Global Change Biology, 23(7), 2768–2782. https://doi.org/10.1111/gcb.13580 Mitsch, W. J., & Gosselink, J. G. (2000). The value of wetlands: Importance of scale and landscape setting. Ecological Economics, 35(1), 25–33. https://doi.org/10.1016/S0921-8009(00)00165-8 84 Oertel, C., Matschullat, J., Zurba, K., Zimmermann, F., & Erasmi, S. (2016). Greenhouse gas emissions from soils—A review. Geochemistry, 76(3), 327–352. https://doi.org/10.1016/j.chemer.2016.04.002 Petrokofsky, G., Kanamaru, H., Achard, F., Goetz, S. J., Joosten, H., Holmgren, P., Lehtonen, A., Menton, M. C., Pullin, A. S., & Wattenbach, M. (2012). Comparison of methods for measuring and assessing carbon stocks and carbon stock changes in terrestrial carbon pools. How do the accuracy and precision of current methods compare? A systematic review protocol. Environmental Evidence, 1(1), 6. https://doi.org/10.1186/2047-2382-1-6 Rey‐Sanchez, C., Arias‐Ortiz, A., Kasak, K., Chu, H., Szutu, D., Verfaillie, J., & Baldocchi, D. (2022). Detecting Hot Spots of Methane Flux Using Footprint‐Weighted Flux Maps. Journal of Geophysical Research: Biogeosciences, 127(8). https://doi.org/10.1029/2022JG006977 Rokni, K., Ahmad, A., Selamat, A., & Hazini, S. (2014). Water Feature Extraction and Change Detection Using Multitemporal Landsat Imagery. Remote Sensing, 6(5), 4173–4189. https://doi.org/10.3390/rs6054173 Roulet, N. T. (2000). Peatlands, carbon storage, greenhouse gases, and the Kyoto Protocol: Prospects and significance for Canada. Wetlands, 20(4), 605–615. https://doi.org/10.1672/0277-5212(2000)020[0605:PCSGGA]2.0.CO;2 Schipper, L. A., & Reddy, K. R. (1994). Methane Production and Emissions from Four Reclaimed and Pristine Wetlands of Southeastern United States. Soil Science Society of America Journal, 58(4), 1270–1275. https://doi.org/10.2136/sssaj1994.03615995005800040039x 85 Sellers, P. J. (1985). Canopy reflectance, photosynthesis and transpiration. International Journal of Remote Sensing, 6(8), 1335–1372. https://doi.org/10.1080/01431168508948283 Singh, K. V., Setia, R., Sahoo, S., Prasad, A., & Pateriya, B. (2015). Evaluation of NDWI and MNDWI for assessment of waterlogging by integrating digital elevation model and groundwater level. Geocarto International, 30(6), 650–661. https://doi.org/10.1080/10106049.2014.965757 Smith, K. A., Ball, T., Conen, F., Dobbie, K. E., Massheder, J., & Rey, A. (2003). Exchange of greenhouse gases between soil and atmosphere: Interactions of soil physical factors and biological processes: Greenhouse gas exchange between soil and atmosphere. European Journal of Soil Science, 54(4), 779–791. https://doi.org/10.1046/j.1351-0754.2003.0567.x Strachan, I. B., Nugent, K. A., Crombie, S., & Bonneville, M.-C. (2015). Carbon dioxide and methane exchange at a cool-temperate freshwater marsh. Environmental Research Letters, 10(6), 065006. https://doi.org/10.1088/1748-9326/10/6/065006 Tanaka, K., & Hashimoto, S. (2006). Plant canopy effects on soil thermal and hydrological properties and soil respiration. Ecological Modelling, 196(1–2), 32–44. https://doi.org/10.1016/j.ecolmodel.2006.01.004 Torres-Alvarado, R., Ramírez-Vives, F., & Fernández, F. J. (2005). Methanogenesis and methane oxidation in wetlands. Implications in the global carbon cycle Metanogénesis y metano-oxidación en humedales. Implicaciones en el ciclo del carbono global. 15(3). Turetsky, M. R., Kotowska, A., Bubier, J., Dise, N. B., Crill, P., Hornibrook, E. R. C., Minkkinen, K., Moore, T. R., Myers-Smith, I. H., Nykänen, H., Olefeldt, D., Rinne, J., Saarnio, S., Shurpali, N., Tuittila, E.-S., Waddington, J. M., White, J. R., Wickland, K. P., & Wilmking, M. (2014). A synthesis of methane emissions from 71 northern, 86 temperate, and subtropical wetlands. Global Change Biology, 20(7), 2183–2197. https://doi.org/10.1111/gcb.12580 Watson, A., Stephen, K. D., Nedwell, D. B., & Arah, J. R. M. (1997). Oxidation of methane in peat: Kinetics of CH4 and O2 removal and the role of plant roots. Soil Biology and Biochemistry, 29(8), 1257–1267. https://doi.org/10.1016/S0038-0717(97)00016-3 Xu, H. (2006). Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. International Journal of Remote Sensing, 27(14), 3025–3033. https://doi.org/10.1080/01431160600589179 87 Appendix A – Meteorological variable gap filling Appendix Figure 1 Pearson correlation tests for MBPPW1 and Hamiota Station. From left to right, 1:1 plots for air temperature, relative humidity, and wind speed along the top row, and soil temperature, incoming shortwave radiation, and wind direction along the bottom row. Values from Hamiota Climate Station are represented along the y-axis, and values from MBPPW1 are represented along the x-axis. The dotted red line represents the 1:1 line, and the orange line is the linear correlation slope between the two data sets. The linear equation and coefficient of correlation R2 are presented in the legend box of each plot. 88 Appendix Figure 2 Pearson correlation tests for MBPPW1 and Oakburn Station. From left to right, 1:1 plots for air temperature, relative humidity, and wind speed along the top row, and soil temperature, incoming shortwave radiation, and wind direction along the bottom row. Values from Oakburn Climate Station are represented along the y-axis, and values from MBPPW1 are represented along the x-axis. The dotted red line represents the 1:1 line, and the orange line is the linear correlation slope between the two data sets. The linear equation and coefficient of correlation R2 are presented in the legend box of each plot. Appendix Figure 1 and Appendix Figure 2 show the 1:1 Pearson correlation analyses between MBPPW1 and Hamiota Climate Station, and MBPPW1 and Oakburn Climate Station respectively. TS stands out as the variable whose linear equation slope and offset, relative to magnitude, deviate farthest from the 1:1 line. This can be attributed to stronger spatial influences from differing canopy and land cover type between sites. Pearson correlation coefficients for TS however are still quite strong at around R2 = 0.921 and 0.899 for Hamiota and Oakburn respectively. The lowest R2 values were found from the two wind variables at 0.771 and 0.787 for wind direction and 0.778 and 0.789 for wind speed However, since wind conditions are 89 expected to be more variable within a regional scale, the lower correlation coefficients for these variables are within an acceptable range. Both climate stations appear to fit the MBPPW1 meteorological data well, so a root mean squared error (RMSE) analysis was performed to compare the goodness of fit from the two climate stations. The results are summarized in Appendix Table 1. Appendix Table 1 Root mean squared error (RMSE) analysis summary. From left to right, RMSE and 𝑅𝑀𝑆𝐸?̅? values are listed for air temperature, relative humidity, wind speed, soil temperature, incoming shortwave radiation, and wind direction. For each meteorological variable, the better performing 𝑅𝑀𝑆𝐸?̅? value between Hamiota and Oakburn are highlighted in green. TA (˚C) RH (%) WS (m/s) TS (˚C) SWIN (W/m2) WD (˚) Hamiota RMSE 1.0499 7.5491 1.2839 3.1121 80.0128 47.0453 𝑹𝑴𝑺𝑬?̅? 0.500 0.105 0.388 0.447 0.510 0.241 Oakburn RMSE 1.589 8.3346 1.7695 3.2545 77.5796 42.5415 𝑹𝑴𝑺𝑬?̅? 0.757 0.116 0.535 0.467 0.494 0.218 RMSE was calculated by the equation, ∑ √(𝑦𝑖−𝑥𝑖)2𝑛𝑛𝑖=1 , where yi is the measured value at MBPPW1 at time i, and xi is the measured value at each respective climate station at time i. RMSE was then divided by the average value of the measurement variable, ?̅?. A smaller RMSE and 𝑅𝑀𝑆𝐸?̅? ratio indicates a better fit. From Table 1, Hamiota’s 𝑅𝑀𝑆𝐸?̅? outperformed Oakburn’s in four out of six tested variables: TA, relative humidity, wind speed, and TS. Oakburn performed notably poorly in TA and wind speed, and only marginally performed better in the remaining incoming shortwave radiation and wind direction variables. From this analysis, it was concluded 90 that Hamiota climate station was more in line with observations at MBPPW1, and was therefore chosen to fill meteorological data gaps at MBPPW1. Only one climate station near MBPPW2, Newdale station, had data available for us to use, so no RMSE goodness of fit comparison was performed for this site. 1:1 correlation tests in Appendix Figure 3 were still used to ensure an acceptable goodness of fit before gap-filling with Newdale measurements. Appendix Figure 3 Pearson correlation tests for MBPPW2 and Newdale Station. From left to right, 1:1 plots for air temperature, relative humidity, and wind speed along the top row, and soil temperature, incoming shortwave radiation, and wind direction along the bottom row. Values from Newdale Climate Station are represented along the y-axis, and values from MBPPW2 are represented along the x-axis. The dotted red line represents the 1:1 line, and the orange line is the linear correlation slope between the two data sets. The linear equation and coefficient of correlation R2 are presented in the legend box of each plot. 91 Pearson correlation coefficients are fairly strong for all observed variables, but a few problems stand out for Newdale. The first is that much like Hamiota and Oakdale, TS and the two wind variables remain weaker fits. The same rationale for why these variables struggled to fit as well as the others at MBPPW1 can also be applied here. The second problem is that the regression slope for incoming shortwave radiation deviated quite strongly from the 1. This deviation was caused from a pyranometer sensor calibration drift at MBPPW2 that was discovered in September 2022, and subsequently corrected in October 2022. Newdale data were concluded to be within the acceptable similarity range and were used to gap-fill missing MBPPW2 meteorological data. 92 Appendix B – Additional figures Appendix Figure 4. Simple linear regression scatterplots for MBPPW1, 2022. 93 Appendix Figure 5. Simple linear regression scatterplots for MBPPW1, 2022. 94 Appendix Table 2 US-Myb r-value data table at α = 0.15 significance level. Pearson correlation coefficient values at the 15% significance level for US-Myb. Positive r-values are highlighted blue, and negative r-values are highlighted orange. US-Myb (α = 0.15) Observation NDVI NDWI NDMI LST 2021-05-08 -0.899 0.837 - -0.873 2021-05-24 - - - -0.852 2021-06-09 -0.753 0.831 - - 2021-06-25 0.780 -0.767 0.952 -0.97 2021-07-11 0.755 -0.726 - -0.944 2021-08-28 0.783 - 0.703 -0.968 2021-09-13 0.762 -0.853 0.865 -0.786 2021-09-29 0.984 -0.987 - -0.793 Appendix Figure 6. Pearson r-value boxplots for US-Myb at the α = 0.15 significance level. 95 Appendix Figure 7 US-Myb (2014) linear regression results. Clustered data from 2014 that were unused due to too many observations yielding p-values far greater than 0.1.
Thesis/Dissertation
eng
Geography
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Knox, Sara
Black, T. Andrew
Characterization and prediction of methane spatial heterogeneity within and beyond a flux tower footprint
Text