SATELLITE TIME-SERIES DATA FOR VEGETATION PHENOLOGY DETECTION AND ENVIRONMENTAL ASSESSMENT IN SOUTHEAST ASIA By Tanita Suepa A DISSERTATION Submitted to Michigan State University in partial fulfillment of requirements for the degree of Geography - Doctor of Philosophy 2013 ABSTRACT SATELLITE TIME-SERIES DATA FOR VEGETATION PHENOLOGY DETECTION AND ENVIRONMENTAL ASSESSMENT IN SOUTHEAST ASIA By Tanita Suepa The relationship between temporal and spatial data is considered the major advantage of remote sensing in research related to biophysical characteristics. With temporally formatted remote sensing products, it is possible to monitor environmental changes as well as global climate change through time and space by analyzing vegetation phenology. Although a number of different methods have been developed to determine the seasonal cycle using time series of vegetation indices, these methods were not designed to explore and monitor changes and trends of vegetation phenology in Southeast Asia (SEA). SEA is adversely affected by impacts of climate change, which causes considerable environmental problems, and the increase in agricultural land conversion and intensification also adds to those problems. Consequently, exploring and monitoring phenological change and environmental impacts are necessary for a better understanding of the ecosystem dynamics and environmental change in this region. This research aimed to investigate inter-annual variability of vegetation phenology and rainfall seasonality, analyze the possible drivers of phenological changes from both climatic and anthropogenic factors, assess the environmental impacts in agricultural areas, and develop an enhanced visualization method for phenological information dissemination. In this research, spatio-temporal patterns of vegetation phenology were analyzed by using MODIS-EVI time series data over the period of 2001-2010. Rainfall seasonality was derived from TRMM daily rainfall rate. Additionally, this research assessed environmental impacts of GHG emissions by using the environmental model (DNDC) to quantify emissions from rice fields in Thailand. Furthermore, a web mapping application was developed to present the output of phenological and environmental analysis with interactive functions. The results revealed that satellite time-series data provided a great opportunity to study regional vegetation variability and internal climatic fluctuation. The EVI and phenological patterns varied spatially according to climate variations and human management. The overall regional mean EVI value in SEA from 2001 to 2010 has gradually decreased and phenological trends appeared to shift towards a later and slightly longer growing season. Regional vegetation dynamics over SEA exhibited patterns associated with major climate events such as El Niño in 2005. The rainy season tended to start early and end late and the length of rainy season was slightly longer. However, the amount of rainfall has decreased from 2001 to 2010. The relationship between phenology and rainfall varied among different ecosystems. Additionally, the local scale results indicated that rainfall is a dominant force of phenological changes in naturally vegetated areas and rainfed croplands, whereas human management is a key factor in heavily agricultural areas with irrigated systems. The results of estimating GHG emissions from rice fields in Thailand demonstrated that human management, climate variation, and physical geography had a significant influence on the change in GHG emissions. In addition, the complexity of spatio-temporal patterns in phenology and related variables were displayed on the visualization system with effective functions and an interactive interface. The information and knowledge in this research are useful for local and regional environmental management and for identifying mitigation strategies in the context of climate change and ecosystem dynamics in this region. ACKNOWLEDGEMENTS The success of my dissertation depended largely on the encouragement, support, and guidance of many important people. I would like to express my deepest appreciation to my advisor, Dr. Jiaguo Qi, who always offers advice and support to his students. He inspired me to do the research in Southeast Asia. Working with him during my graduate study provides more valuable experiences and it helped me to develop my research ideas and intellectual skills. He always provided very useful comments, remarks, and engagement through the learning process of this dissertation. Without his guidance and persistent help, this dissertation would not have been possible. I give a special gratitude to my committee members, Dr. Joseph Messina, Dr. Kirk Goldsberry, and Dr. Sasha Kravchenko. Their contribution in sharing invaluable knowledge, wisdom, advice, and comments, as well as precious time, with me, helped me go through the dissertation process. Dr. Joe brought the technical and methodological perspectives to my research and I would like to thank for his professional advice. I wish to thank Dr. Sasha for her kindness and her expertise in statistical comments. I am thankful for Dr. Kirk’s comments and suggestions on this research. Although he has lived in Boston for the last two years, he helped me to shape and improve my work from the very beginning of my study and he inspires me to learn and understand about cartography, which is very useful for my career. I would also like to acknowledge with much appreciation to my scholarship sponsors, the Royal Thai Government, the Ministry of Science and Technology, and the Geo-informatics and Space Technology Development Agency (GISTDA) for their financial support, which opened my bright vision on remote sensing and cartography and made me get this prestigious degree. iv The funding from the Department of Geography, the Graduate School, and the Office for International Student and Scholars at Michigan State University, as well as Kenneth E. and Marie J. Corey Research Enrichment Fund Awards, provided the opportunity for me to complete this research. Collecting field survey, getting valuable data, and attending international conferences were possible with these funding. I would like to extend my gratitude to the following people and organizations, Dr. Arika Bridhikitti (Mahasarakham University), Dr. Kruamas Smakgahn (Kasetsart University), the Office of Agricultural Economics (the Ministry of Agriculture and Cooperatives), the Thai Meteorological Department, for supporting the data and information for my research. Furthermore, I sincerely thank my friends in Thailand, especially my friends at GISTDA, for supporting me and encouraging me as well as helping me to get the data from Thai organizations. Without them, it might have taken more time to get the data to complete this research. A special thanks goes to Dr. Feng Zhang, Lanzhou University, who provided very helpful comments and technical advice as well as guided me through difficulties in the DNDC model. I offer my heartfelt thanks to my friends, Jenni Gronseth, Alisa Pamela Lertpratchya, and Chanaichon Damsri, who helped me to edit my dissertation and improve my writing skills. I sincerely appreciate all they have done for me. I would like to express my gratitude to the Department of Geography and the Center for Global Change and Earth Observations (CGCEO) at Michigan State University for teaching me the knowledge necessary to be a good geographer and providing the wonderful academic facility. I would like to show my greatest appreciation to Sharon Ruggles, Graduate Secretary of Department of Geography, for everything that she has done for me over the past five years. I must thank her for being so kind and active to respond my questions and problems. Many thanks v go to Jean Lyle Lepard, Kathleen Mills, Cameron Williams for their assistance during my fives year at Michigan State University. I am also grateful to all my friends at Remote Sensing Lab and my Thai friends in the US. who made it enjoyable to do my research and helped me recover from homesickness. A special thanks goes out to my brother and my friend, Siam Lawawirojwong, for helping me to solve many technical issues, especially Python programming. Without him, my data processing might have taken more than a year to complete. Thank you for standing by me, supporting me, and keeping my spirits up through the stressful times. I could not have completed this research without his support. Finally, my deepest thank to my family, my mother, father, and sister, for everything they have done to help me reach my goals. They have supported me and encouraged me throughout entire process. Their unending love and support made this achievement possible. I will be grateful forever for their love. vi TABLE OF CONTENTS LIST OF TABLES………………………...…..………………………………………….. xi LIST OF FIGURES………………………...…………………………………………….. xii ABBREVIATIONS………………………...…..…………………………………………. xvi CHAPTER 1 INTRODUCTION…………………………………….……………...….... 1.1 Research Problems…………………………....…..………………....……………… 1.2 Conceptual Framework…………………………………………...………………… 1.3 Research Objectives and Research Questions………………………………………. 1.3.1 Objectives…………………………………………...……………………… 1.3.2 Research Questions…………………………………………...…………….. 1 1 6 8 8 8 CHAPTER 2 LITERATURE REVIEW………………………………….……………... 2.1 Satellite Time-series Data for Phenological Information.………………………….. 2.1.1 Satellite Time-series Data and Vegetation Phenology………….………….. 2.1.2 Remote Sensing Approaches…………………………………….…………. 2.1.3 Vegetation Index…………………………………………………….……... 2.1.4 Phenological Information Derived from Satellite Time-series Data….……. 1) Growing Season Profile and Phenological Parameters………………….. 2) Data Smoothing and Filtering…………………………………………… 3) Extraction of Phenological Parameters………………………………….. 4) Trend Analysis for Satellite Time-series Data…………………………... 2.2 Phenology and Climate Variability…………………………………………………. 2.2.1 Relationship between Phenology and Climate…….……………….………. 2.2.2 Satellite Time-series Data for Vegetation Phenology Response to Rainfall Change……………………………………………………………………… 2.3 Environmental Model: The DNDC Model………………………………………….. 2.3.1 Global Climate Change and Greenhouse Gases Emissions………….…….. 2.3.2 DNDC Model and its Capacities……………………………………….…... 2.3.3 GHG Emissions in SEA and Thailand………………………………….….. 2.4 Visualization for Satellite Phenological Information……………………………….. 2.4.1 Remote Sensing Time-series Data and Visualization……………….……... 2.4.2 Visual Exploration……………………………………………..…………… 2.4.3 Cartographic Design for Satellite Time-series Data………….……………. 1) Internet Atlas…………………………………………………….………. 2) Thematic Map…………………………………………………….……... 3) Bivariate Map…………………………………………………….……… 4) Map Animation…………………………………………………….……. 5) Web Mapping Application: A Tool to Disseminate Geographic Data….. 10 10 10 14 15 18 18 19 24 26 28 28 vii 30 33 33 35 37 44 44 45 46 46 47 48 49 51 CHAPTER 3 RESEARCH METHODOLOGY………………..………….……………. 3.1 Study Area…………………………………………………………….…………….. 3.2 Data Sources………………………………………………………………………… 3.3 Data Pre-processing…………………………………………………………………. 3.4 Fitted Function and Phenological Extraction……………………………………….. 3.4.1 Data Filtering……………………………….………………………………. 3.4.2 Phenological Parameters………………………….………………………... 3.5 Rainfall Seasonality Extraction……………………………………………………... 3.6 Spatial Variation and Trend Analysis………………………………………………. 3.6.1 Trend Analysis……………………………….…………………………….. 3.6.2 Spatio-temporal Patterns, Changes and Trends of EVI Profile…….………. 3.6.3 Spatio-temporal Patterns, Changes and Trends of Phenological Parameters. 3.6.4 Spatio-temporal Patterns, Changes and Trends of Rainfall Seasonality….... 3.6.5 Relationship between Phenology and Rainfall Seasonality……….….......... 3.6.6 Analyze the Drivers of Phenology Changes…………………….…………. 3.6.7 Comparison between MODIS Phenology and Field Observations and TRMM Rainfall and Station Rainfall………………………………………. 3.7 Modeling DNDC……………………………………………………………………. 3.7.1 Model Description………………………………………………….………. 3.7.2 Study Area……………………………………………………….…………. 3.7.3 Database Development and Input Parameters………………….…………... 3.7.4 Spatio-temporal Patterns and Changes of GHG Emissions under Different Scenarios………………………………………………………….………… 3.7.5 Comparison of DNDC Results to IPCC Approach and Thailand Research…………………………………………………………………….. 3.8 The Interactive Phenological Atlas for SEA (IPA)………………………………..... 3.8.1 Conceptual Framework…………………………………………………….. 3.8.2 Cartographic Design for Satellite Time-series Data……………………….. 3.8.3 Visualization Implementation……………………………………………… 54 54 57 58 60 61 64 66 67 67 70 70 70 71 71 CHAPTER 4 SPATIO-TEMPORAL VARIATION OF VEGETATION PHENOLOGY IN SEA…………………………………………………… 4.1 Introduction………………………………………………………………………… 4.2 Seasonal Characteristics of Vegetation Profiles and Land Use/ Land Cover Types in SEA………………………………………………………………………………. 4.3 Spatio-temporal Variation of EVI………………………………………………….. 4.3.1 Spatial Patterns of Interannual EVI………………………………………… 4.3.2 Trend Analysis……………………………………………………………… 4.4 Spatio-temporal Variation of Phenological Parameters……………………………. 4.4.1 Spatial Patterns of Interannual Phenology………………………………….. - The First Growing Season (Season 1) - The Second Growing Season (Season 2) 4.4.2 Trend Analysis……………………………………………………………… - The First Growing Season (Season 1) …………………………………. - The Second Growing Season…………………………………………… 4.4.3 Trend from Bivariate Mapping……………………………………………... 98 viii 72 74 74 78 78 86 88 88 88 91 93 98 99 107 107 108 110 110 114 117 118 118 121 133 - The Start and End of the Growing Season……………………………... - The Amplitude and Length of the Growing Season……………………. - The Large Integral and Length of the Growing Season………………... 4.5 Comparison between MODIS Phenology and Field Observations………………… 4.6 Discussion and Conclusions………………………………………………………... 4.6.1 Spatio-temporal Variation of Vegetation Phenology in SEA…………….... 4.6.2 Challenges and Limitations………………………………………………… - Data Quality…………………………………………………….………. - Validation…………………………………………………….………… - Resolution…………………………………………………….………… 134 135 135 140 142 142 144 145 147 148 CHAPTER 5 UNDERSTANDING THE DRIVING FORCES OF PHENOLOGICAL CHANGES.…………………………………………. 5.1 Introduction………………………………………………………………………… 5.2 Response of Seasonal Vegetation Dynamics to Climate Variations in SEA……….………………………………………………………………………… 5.2.1 Spatial-temporal Variation on Rainfall Seasonality………………………... 5.2.2 Comparison between Mean Annual EVI and Mean Annual Rainfall……… 5.2.3 Relationship between Seasonal Rainfall Fluctuations and Phenological Parameters………………………………………………………………….. 5.2.4 Rate of Phenological Changes with Respect to Rainfall Seasonality………………………………………………………………...... 5.3 Driver of Phenological Change: Local Scale Analysis…………………………….. 5.3.1 Site I and II in Northern Thailand………………………………………….. 5.3.2 Site III and IV in Central Thailand…………………………………………. 5.4 Comparison between TRMM Rainfall and Station Rainfall……………………….. 5.5 Discussion and Conclusions………………………………………………………... 5.5.1 Spatio-temporal Variation of Vegetation Phenology and Rainfall Seasonality ………………………………………………………………… 5.5.2 The Influence of Climate Change on Phenological Patterns……………….. 5.5.3 The Influence of Human Management on Phenological Patterns………….. 150 150 151 151 152 152 153 159 160 162 169 172 172 173 177 CHAPTER 6 ENVIRONMENTAL IMPACTS WITH BIOGEOCHEMICAL MODEL: RICE FIELDS IN THAILAND………………………………. 179 6.1 Introduction………………………………………………………………………… 179 6.2 Spatio-temporal Patterns of GHG Emissions under Different Scenarios…………... 180 6.2.1 Scenario 1: Actual Input Data for 2002 and 2010………………………….. 180 6.2.2 Scenario 2: Fertilizer Effect 2002 and 2010……………..…………………. 183 6.2.3 Scenario 3: Climate Effect…………………………..……………………... 191 6.2.4 Scenario 4: Phenology Effect………………………..……………………... 194 6.2.5 Scenario 5: Climate and Phenology Effect………………………..……….. 195 6.2.6 Global Warming Potential of Lopburi Province………………………..….. 198 6.3 Comparison to IPCC Estimation and Research in Thailand……………………….. 199 6.4 Discussion and Conclusions………………………………………………………... 202 6.4.1 Effect of Soil Properties on GHG Emissions……………………………..... 203 6.4.2 Effect of Management Practices on GHG Emissions……………………… 204 ix 6.4.3 6.4.4 6.4.5 Effect of Climate and Phenology on GHG Emissions……………………... Challenges and Limitations……………………………................................ Potential of DNDC Model…………………………….................................. 207 208 209 CHAPTER 7 THE INTERACTIVE PHENOLOGICAL ATLAS FOR SEA………… 7.1 Introduction………………………………………………………………………… 7.2 Design and Implementation of IPA………………………………………………… 7.2.1 Graphical User Interface (GUI) Design…………………………………….. 7.2.2 Map Design…………………………………………………………………. 7.2.3 Creating Map Service………………………………………………………. 7.2.4 Creating Map Animation…………………………………………………… 7.2.5 Establishing and Customizing the Prototype of IPA……………………….. 7.3 Capabilities of IPA……………………………………………………………......... 7.3.1 Accessibility………………………………………………………………... 7.3.2 Data…………………………………………………………………………. 7.3.3 Functionality………………………………………………………………... 1) Basic Functions…………………………………………………….……. 2) Advance Functions…………………………………………………….… 7.4 Discussion and Conclusions………………………………………………………... 7.4.1 Achievement: The Concept of Visual Exploration in IPA………………..... 7.4.2 Insights: Challenges and Limitations. ..………………................................. 213 213 214 214 218 221 221 222 223 223 223 223 224 224 226 226 227 CHAPTER 8 CONCLUSIONS AND FUTURE ENVISIONS.………………………… 8.1 Major Findings……………………………………………………………………… 8.1.1 Spatio-temporal Variation of Vegetation Phenology and Rainfall Seasonality………………………………………………………………….. 8.1.2 Spatio-temporal patterns of GHG Emissions from the DNDC Model under Different Scenarios…………………………………………………………. 8.1.3 The Interactive Phenological Atlas for SEA……………………………….. 8.2 Challenges and Limitations…………………………………………………………. 8.3 Future Research……………………………………………………………………... 230 230 APPENDICES…………………………………………………………………………… Appendix A: EVI and Phenological Parameters 2001-2010……………………………….. Appendix B: Rainfall Seasonality 2001-2010……………………………………………… Appendix C: Data Quality Assessment Map……………………………………………….. Appendix D: Multi-year Simulation of GHG Emissions from DNDC Model (Scenario 3-5)……………………………..………………………………….. 241 242 246 250 REFERENCES……………………………………………………………………………. 261 x 230 233 235 237 238 252 LIST OF TABLES Table 3.1 Summary of data requirements in this research…………………………………. 57 Table 3.2 TIMESAT input parameters……………………………………………………... 62 Table 3.3 Phenological parameters………………………………………………………… 65 Table 3.4 Input parameters for DNDC model……………………………………………… 84 Table 3.5 Five simulated GHG emissions from rice field in Thailand under different scenarios…………………………………………………………………………. 87 Table 3.6 Data and map contents in visualization system………………………………….. 90 Table 5.1 Land use change for site I Sukhothai province in 2002 and 2009………………. 161 Table 5.2 Land use change for site II Uttaradit province in 2001 and 2009……………….. 162 Table 5.3 Land use change for site III Suphanburi province in 2000 and 2010…………… 164 Table 5.4 Land use change for site IV Lopburi province in 2002 and 2010……………….. 164 Table 6.1 Emission rate in scenario 1……………………………………………………… 182 Table 6.2 Total emission in scenario 1……………………………………………………... 183 Table 6.3 Emission rate in scenario 2……………………………………………………… 190 Table 6.4 Total emission in scenario 2……………………………………………………... 190 Table 6.5 Comparison between modeled estimation and IPCC approach…………………. 201 Table 6.6 Comparison between modeled estimation and Thailand standard scaling factors. 202 Table 6.7 Comparison between modeled estimation and standard emission rate in Lopburi province………………………………………………………………………….. 202 xi LIST OF FIGURES Figure 1.1 Conceptual framework………………………………………………………….. 6 Figure 2.1 A simple NDVI profile for a typical patch of vegetation………………………. 20 Figure 2.2 The results of data smoothing applied to EVI time-series, the start and end of the seasons are marked with filled circles……………………………………… 20 Figure 3.1 Research methodology………………………………………………………….. 55 Figure 3.2 Study area, the Indochina Peninsula in SEA…………………………………… 56 Figure 3.3 Flowchart of phenology and rainfall analysis…………………………………... 59 Figure 3.4 Three year time-series data to extract phonological parameters from one year... 61 Figure 3.5 Testing the range of growing season in Thailand………………………………. 66 Figure 3.6 Field points in central and northeastern Thailand…………………………………... 73 Figure 3.7 DNDC model…………………………………………………………………… 76 Figure 3.8 DNDC flow chart……………………………………………………………….. 77 Figure 3.9a Rice cropping frequency in Lopburi province 2002…………………………... 79 Figure 3.9b Rice cropping frequency in Lopburi province 2010…………………………... 80 Figure 3.10 Interactive phenological atlas for SEA………………………………………... 89 Figure 3.11 Web mapping architecture…………………………………………………….. 94 Figure 4.1 Land cover in SEA from GLOBCOVER in 2005……………………………… 100 Figure 4.2 Rice cropping system in SEA 2007…………………………………………….. 102 Figure 4.3 Rice Cropping Frequency in SEA 2007………………………………………... 103 Figure 4.4 Seasonal EVI profile patterns…………………………………………………... 105 Figure 4.5 Seasonal EVI profile changes…………………………………………………... 106 xii Figure 4.6 (a) 10-year mean EVI (2001-2010) and (b) EVI trend (2001-2010)…………… 111 Figure 4.7 EVI mean deviation 2001-2010………………………………………………… 111 Figure 4.8 EVI trend of hotspots in forest areas…………………………………………… 112 Figure 4.9 EVI trend of hotspots in agricultural areas……………………………………... 112 Figure 4.10 10-year mean phenology, the first growing season (2001-2010)……………... 113 Figure 4.11 10-year mean phenology, the second growing season (2001-2010)…............... 122 Figure 4.12 Phenology trend, the first growing season (2001-2010)………………………. 124 Figure 4.13 Phenology trend of hotspots in forest areas…………………………………… 126 Figure 4.14a Phenology trend of hotspots in agricultural areas (the first growing season)... 127 Figure 4.14b Phenology trend of hotspots in agricultural areas (the second growing season)………………………………………………………………………... 129 Figure 4.15 Phenology trends, the second growing season (2001-2010)…………………... 131 Figure 4.16 Bivariate map: mean different between the start and the end dates (2001-2010)……………………………………………………………………. 137 Figure 4.17 Bivariate map: mean different between the length and the amplitude (2001-2010)……………………………………………………………………. 138 Figure 4.18 Bivariate map: mean different between the length and the large integral (2001-2010)……………………………………………………………………. 139 Figure 4.19 Comparison between MODIS phenology and field data in Thailand…………. 141 Figure 5.1 10-year mean rainfall (2001-2010)……………………….…………………….. 154 Figure 5.2 Rainfall trends (2001-2010)……………………….............................................. 155 Figure 5.3 Mean annual EVI and mean annual rainfall…………………………………..... 156 Figure 5.4 Relationship between phenology and rainfall seasonality……………………… 157 Figure 5.5 Rate of phenological change with respect to rainfall seasonality……………..... 158 Figure 5.6 Land use changes in site I and II in northern Thailand………………………..... 165 xiii Figure 5.7 Land use changes in site III and IV in central Thailand………………………... 167 Figure 5.8 Comparison annual rainfall between TRMM rainfall and station rainfall in Thailand……………………………………………………………………….. 170 Figure 5.9 Comparison rainfall parameters between TRMM rainfall and station rainfall in Thailand………………………………………………………………………… 171 Figure 6.1 GHG emission rate in scenario 1……………………………………………….. 184 Figure 6.2 Total GHG emission in scenario 1……………………………………………… 186 Figure 6.3 Comparison GHG emission rate between scenario 1 (S1) and scenario 2 (S2)... 188 Figure 6.4 Soil texture and % clay in Lopburi province…………………………………… 189 Figure 6.5 GHG emission rate in scenario 2……………………………………………….. 192 Figure 6.6 Comparison between climate input data and GHG emission rate in scenario 3... 196 Figure 6.7 Comparison between phenology input data and GHG emission rate in scenario 4……………………………………………………………………….. 197 Figure 6.8 GHG emission rate in scenario 5……………………………………………….. 198 Figure 7.1 IPA interface and components………………………………………………….. 216 Figure A1 EVI and phenological parameters 2001-2010…………………………………... 242 Figure A2 Start date of the Growing Season 2001-2010 (first growing season) ………….. 243 Figure A3 End date of the Growing Season 2001-2010 (first growing season) …………... 244 Figure A4 Length of the Growing Season 2001-2010 (first growing season) …………….. 245 Figure B1 Mean Annual Rainfall 2001-2010………………………………………………. 246 Figure B2 Start Date of the Rainy Season 2001-2010……………………………………... 247 Figure B3 End Date of the Rainy Season 2001-2010……………………………………… 248 Figure B4 Length of the Rainy Season 2001-2010………………………………………… 249 Figure C1 MODIS EVI data quality Assessment 2001-2010……………………………… 250 xiv Figure C2 MODIS EVI data quality assessment (10-year mean) …………………………. 251 Figure D1 CH4 emission rate in scenario 3………………………………………………... 252 Figure D2 CO2 emission rate in scenario 3………………………………………………... 253 Figure D3 N2O emission rate in scenario 3………………………………………………... 254 Figure D4 CH4 emission rate in scenario 4………………………………………………... 255 Figure D5 CO2 emission rate in scenario 4………………………………………………... 256 Figure D6 N2O emission rate in scenario 4………………………………………………... 257 Figure D7 CH4 emission rate in scenario 5………………………………………………... 258 Figure D8 CO2 emission rate in scenario 5………………………………………………... 259 Figure D9 N2O emission rate in scenario 5………………………………………………... 260 xv ABBREVIATIONS ADB Asian Development Bank AG Asymmetric Gaussian Functions AP Ammonium phosphate API Application Programming Interface AVHRR Advanced Very High Resolution Radio Meter from the National Oceanic and Atmospheric Administration’s satellite C Carbon CH4 Methane CO2 Carbon dioxide DAP days after planting DBH days before harvest DBP days before planting DL The Double Logistic Function DNDC DeNitrification and DeComposition model DOY The day of the year of rainy season EGS End of the growing season Eh Soil Redox Potential ENSO El Niño-Southern Oscillation events EPA United States Environmental Protection Agency ERS The end of rainy season EVI The Enhanced Vegetation Index xvi FAO Food and Agriculture Organization of the United Nations FT Fourier Transform GHG Greenhouse gases emissions GIS Geographic Information System GUI Graphical User Interface GWP Global Warming Potential IPA Interactive Phenological Atlas for SEA IPCC Intergovernment Panel on Climate Change LGS Length of growing season LP DAAC Land Processes Distributed Active Archive Center LRS The length of rainy season LSP Land surface phenology MK Mann-Kendall MODIS The Moderate Resolution Imaging Spectroradiometer mV Millivolt N Nitrogen NH3 Ammonia N2O Nitrous oxide NDVI The Normalized Difference Vegetation Index NIC National Intelligence Council NIR Near Infrared NO Nitrogen oxide NOAA National Oceanic and Atmospheric Administration xvii NOx oxides of nitrogen PN Phenology R coefficient of correlation RD relative deviation SEA Southeast Asia SG Savitzky-Golay Filtering SGS Start of the growing season SK Seasonal Kendall SOC Soil organic carbon SRS The start of rainy season TMD Thailand Meteorological Department TRMM Tropical Rainfall Measuring Mission UMN University of Minnesota VGT the Vegetation Instrument Sensors of SPOT VI Vegetation Index WMS Web map services xviii CHAPTER 1 INTRODUCTION 1.1 Research Problems The relationship between temporal and spatial data is considered the major advantage of remote sensing in research related to biophysical characteristics. Due to the strong interactions between global climate change and vegetation phenology, research in plant phenology has brought on intense interest in the context of climate change. Remote sensing data can contribute significantly to this research. With temporally formatted remote sensing products, it is possible to monitor environmental changes through time and space by analyzing vegetation processes. Additionally, remote sensing technology has changed the observation of plant phenology from points (phenological observation stations) to coverages (macro regions) for better phenology observation across regions, countries, continents, and even across the globe (Zhang et al., 2006). Vegetation phenology is the study of recurring patterns of vegetation growth and development, as well as their connection to climate (White et al., 1997). Phenological properties, such as the timing and rate of green-up, amplitude and duration of the growing season, and timing and rate of senescence of plant classes, have become an emerging indicator of global environmental changes. Land surface phenology (LSP) can be observed from remote sensing to monitor the seasonal pattern of the spatio-temporal variation in the vegetated land surfaces (White & Nemani, 2006; de Beurs & Henebry, 2010). Therefore, LSP is a key indicator of ecosystem dynamics under a changing environment (Xiao et al., 2009). Plant phenology is sensitive to changes in weather and climate. Changes in plant phenology can affect the carbon cycle, water cycle, and energy fluxes through photosynthesis and evapotranspiration (Xiao et al., 2009), which consequently may influence food security, 1 water resources availability, and climate. The length or magnitude of the plant growing season may change only slightly, but this could result in large changes in annual gross primary production. Furthermore, shifts in the timing of plant activity (e.g., the start of the growing season) provide evidence that species and ecosystems are being influenced by global environmental change (Reed et al., 2009). Accordingly, it is useful to characterize land cover phenology and understand how phenology responds to interannual variability, climate, and land use. Multi-temporal remote sensing data provide opportunity to characterize LSP at the regional level (Reed et al., 2003). In contrast to AVHRR (Advanced Very High Resolution Radio Meter from the National Oceanic and Atmospheric Administration’s satellite) data as well as SPOT (Système Pour l'Observation de la Terre) and LANDSAT, MODIS (The Moderate Resolution Imaging Spectroradiometer) provides an improved basis for monitoring global ecosystem dynamics with appropriate spatial and temporal resolutions and substantially improved geometric and radiometric properties (Zhang et al., 2006). Additionally, The MODIS EVI ( Enhanced Vegetation Index) is designed to enhance the vegetation signal with improved sensitivity in high biomass regions and reduce atmospheric effects and soil background as well as reduce the smoke of biomass burning in tropical areas (Huete et al., 2002). Therefore, EVI is appropriate for studying phenological dynamics in tropical zones with several growing seasons during the year. Although a number of different methods have been developed to determine the seasonal cycle of vegetation using time-series of vegetation indices, these studies have not explored and monitored changes and trends of vegetation phenology in Southeast Asia (SEA), which is adversely affected by impacts of climate change. SEA is one of the world’s regions that is most 2 vulnerable to the impacts of climate change because of its unique economic and social characteristics, extensive coastal regions, mostly tropical climate, and growing population as well as high dependence upon natural resources for development (IPCC, 2007). As one of the world’s most dynamic regions, the rapid economic growth in the past few decades has been dependent on agriculture. The effects of climate change due to heat stress, water stress, extreme climate events, and a rise in sea levels have caused the agriculture in SEA to be under considerable environmental pressure (IPCC, 2007; ADB, 2009). In the meantime, land use change and rapid economic development have exacerbated the effects of climate change on human and natural systems in this region. Subsequently, it is of great interest to assess the environmental change of this region to better understand the environmental consequences of climate change and of human induced land use land cover change. Changes and trends of vegetation phenology can provide strong scientific evidence for ecosystem dynamics. Furthermore, the relationship between phenology and climate variability is necessary to monitor and assess environmental changes in this region. This relationship is a key to understanding the influence of climate change on biophysical characteristics. In tropical regions, ecosystems are less sensitive to temperature but are dependent on rainfall to trigger the emergence of green leaves and control vegetation growth duration (Kramer et al., 2000; Cleland et al., 2007). More importantly, extreme climate events, i.e., floods and droughts, also influence phenology in these regions (Zhang et al., 2005). Therefore, it is necessary to examine the response of phenology to rainfall variation in SEA in order to understand the sensitivity of various vegetation formations to climate variability and environmental change in this region. 3 Additionally, SEA is considered to be a major source of food production in the world, particularly of rice that is mostly grown by Vietnam and Thailand (FAO, 2007). Increasing demand for food and industrial crops has led to increases in agricultural land conversion and intensification, generating considerable environmental problems. This land use change has direct environmental impacts, particularly the emission of greenhouse gases (GHG). In comparison to global GHG emissions, SEA contributed 12% of the world’s GHG emissions in 2000 and showed an increase of 27% from 1990, which was greater than the global average (ADB, 2009). In addition, SEA countries employ advanced farming practices (e.g., irrigation, tillage, fertilizers) to increase yields and agricultural production. These processes have significant impacts on the environment resulting in increasing GHG emissions (Carbon dioxide (CO2), methane (CH4), and nitrous oxide (N2O)), particularly from rice fields. In order to assess the impact of environmental change, particularly in agricultural areas which is the major source of emissions in this region, an environmental model is needed (ADB, 2009). Although there are several models for estimating GHG emissions, the DNDC (DeNitrification and DeComposition) model has proven to be capable of quantifying emissions in SEA. This is because DNDC has been continuously modified and calibrated as well as validated in Thailand (Cai et al., 2003; Smakgahn et al., 2009). This model can be applied for a wide variety of crop types in both upland and wetland agricultural ecosystems. In addition, DNDC includes CH4 emissions and considers water management in rice cultivation. In addition to the analyses of spatio-temporal patterns, phenological changes, and environmental assessments, a visualization system is essential to provide the understanding and knowledge of how these changes are linked to environmental consequences. Due to the vast 4 amount of information that results from phenological and environmental analyses, a visualization system (e.g., map animation) has become indispensable to represent these spatio-temporal processes and to simulate the dynamics of environmental systems. The visualization system with a dynamic display and interactive interface enables users not only to detect changes but also to understand meanings encoded in transitions within these dynamic displays (Fish et al., 2011). Although phenological information derived from satellite time-series data have provided a great deal of knowledge in environmental change, most research has focused on the northern high latitudes such as the US and China. Additionally, such methods do not account for specific ecosystem characteristics, such as multiple growth cycles and the tropical forest, and they have not addressed the phenology changes and trends in SEA. Therefore, SEA is an ideal case study for exploring and monitoring patterns and trends of vegetation phenological changes in order to indicate the effects of climate variability as well as anthropogenic influences in this region. According to these important issues in SEA, this research developed effective methods to study phenological characteristics in SEA. This study provided information on the spatial distribution and the temporal trends of phenological changes, as well as analyzed the possible drivers of changes between anthropogenically driven land cover change and interannual climatic fluctuations. The results indicate interannual variations of phenology, which correspond to land cover and climate variability. The spatial distributions and changes of GHG emissions were also investigated by using an environmental model. Additionally, a visualization system was developed to explore and monitor patterns and changes of vegetation phenology in this region. The relationship between phenological changes and climate variability, as well as environmental impacts, could enhance the understanding of environmental consequences and ecological changes in SEA. A better understanding of regional phenology and its responses to climate as 5 well as environmental impacts is important for global environmental impact and ecosystem dynamics. This information and knowledge can further identify mitigation strategies and decision-making in the context of environment and climate change for this region. 1.2 Conceptual Framework This research explored ecosystem dynamics in SEA through LSP, the relationship between LSP and climate data (rainfall), environmental impacts (CO2, CH4, and N2O), and the visualization system (Figure 1.1). Figure 1.1 Conceptual framework “For interpretation of the refereces to color in this and all other figures, the reader is referred to the electronic version of this dissertation” LSP is the seasonal pattern of the spatio-temporal variation in the vegetated land surfaces observed from remote sensing (White & Nemani, 2006; de Beurs & Henebry, 2010). The temporal and spatial characteristics of remote sensing make it possible to study vegetation conditions and phenological changes through time and space. Generally, a vegetation index from satellite time-series data is used to explore LSP. In this study, EVI was applied to track 6 vegetation phenology to determine spatio-temporal patterns as well as trends and changes in SEA. Additionally, this research analyzed the causes or drivers of phenological changes by addressing the relationship between phenology and climate variability. This relationship is a key to understanding the influence of climate change on biophysical characteristics. Differences in phenological trends between ecosystems and years are suggestive of ecosystem susceptibility to future climate change. This research used rainfall data as the key factor to explore phenology in response to climate. Land use data were also analyzed to indicate the possible causes of change in phenology by human management. Changes in vegetation phenology and land use have an important impact on environment, particularly change in GHG emissions. To assess the environmental impacts, the DNDC environmental model was applied to quantify GHG emissions (CO2, CH4, and N2O) from rice fields in Thailand. The last concept is the visualization system that is essential to provide the understanding and knowledge of phenological changes and environmental assessment. The visualization system in this research (Interactive Phenological Atlas for SEA) uses a web mapping application to represent the outputs of phenological and environmental analysis in an Internet atlas with interactive functions and powerful techniques such as map animation and bivariate mapping. This system will enable users to increase their understanding of temporal and spatial changes of geographic phenomena and facilitate environmental monitoring and management tasks. These concepts—land surface phenology, the relationship between LSP and climate data, environmental impacts, and the visualization system—were integrated to understand ecosystem 7 dynamics and environmental consequences in SEA. This knowledge will be useful for decisionmaking and environmental management in the near-future. 1.3 Research Objectives and Research Questions 1.3.1 Objectives 1. Determine the spatial characteristics of the vegetation phenology in SEA by extracting seasonal parameters related to the growing seasons; 2. Identify seasonal and inter-annual variability of vegetation phenology and analyze the distribution and the patterns of vegetation phenological changes; 3. Analyze the causes/drivers of phenological changes from both climatic and anthropogenic factors 4. Assess the environmental impacts due to phenological change by using DNDC models for SEA; 5. Develop an improved visualization method for phenological information dissemination. 1.3.2 Research Questions Objective 1: Determine the spatial characteristics of the vegetation phenology in SEA by extracting seasonal parameters related to the growing seasons 1. What are the spatio-temporal patterns of vegetation phenology in SEA, particularly in multiple growth cycle areas? Objective 2: Identify seasonal and inter-annual variability of vegetation phenology and analyze the distribution and the patterns of vegetation phenological changes 1. Does/how does vegetation phenological distribution change over time and across the space at regional and local scales? 8 2. What are the trends of vegetation phenological change in this region? 3. Do the changes of vegetation phenology have magnitudes of difference between agricultural areas and naturally vegetated areas? Objective 3: Analyze the causes/drivers of changes both climatic and anthropogenic factors 1. What are the spatio-temporal patterns and trends of rainfall seasonality? 2. Do the changes and trends of vegetation phenology in SEA have a significant correlation with rainfall change? 3. Are the changes and trends of vegetation phenology related to land use/land cover variation? Objective 4: Assess the environmental impacts due to phenological change by using DNDC models for SEA 1. What are the spatial distributions and changes of GHG emissions (CO2, N2O, CH4) in agricultural areas? 2. Are the changes and distributions of GHG emissions correlated with phenological changes and climate variation? Objective 5: Develop an improved visualization method for phenological information dissemination 1. What is the appropriate technique to compile and organize all phenological information into an effective visualization system? 2. What visualization tools reveal geographical distribution and spatio-temporal patterns of the phenological changes? 9 CHAPTER 2 LITERATURE REVIEW 2.1 Satellite Time-series Data for Phenological Information 2.1.1 Satellite Time-series Data and Vegetation Phenology Temporally formatted remote sensing products enable the analysis of vegetation time series to monitor the process and the pattern of environmental changes through time and space. The relationship between the temporal and the spatial data is considered the major advantage of remote sensing data. These data with various spatial resolution and high temporal resolution provide effectively monitoring large-scale geographic phenomenon, particularly the analysis of the vegetation process. Phenology is a key indicator of ecosystem dynamics under a changing environment (Xiao et al., 2009). According to Schwartz (2003), plant phenology is the study of periodic biological events in the plant world or the timing of different stages of the vegetation seasonal cycle, such as leaf unfolding, first bloom, and leaf fall, which are influenced by the environment and climate. Phenological information—such as the timing and rate of green-up, amplitude and duration of the growing season, and timing and rate of senescence of plant classes—is a key to understand climate variability and trends on vegetation which are important factors in the global change sciences (Hermance, 2007). Land surface phenology (LSP) is defined as the seasonal pattern of the spatiotemporal variation in the vegetated land surfaces observed from remote sensing (White & Nemani, 2006; de Beurs & Henebry, 2010). This definition refers to aggregated information of plant phenology at the spatial resolution of satellite sensors. Therefore, LSP is different from the traditional definition, which refers to specific life cycle events using in situ observations of 10 individual plants or species (Tan et al., 2011). LSP is usually addressed through vegetation indices of satellite time-series data and it can present the advantage of a global coverage. Xiao et al. (2009) mentioned that LSP is sensitive to changes in climate. The carbon cycle, water cycle and energy fluxes change according to changes in phenology of plant through photosynthesis and evapotranspiration. Therefore, these effects have an influence on food security, water resources availability and climate. Phenology can change slightly, such as the length or magnitude of the plant growing season, could result in large changes in annual gross primary production. Furthermore, shifts in the timing of plant activity provide evidence that species and ecosystems are being influenced by global environmental change (Reed et al., 2009). Differences in phenology trends between ecosystems can be used to assess responsiveness to climate change. The degree of interannual variability of phenology, particularly during severe dry and wet years, is suggestive of ecosystem susceptibility to future climate change (Bradley & Mustard, 2008). Accordingly, it is important to characterize land cover phenology and understand how phenology responds to interannual variability, climate, and land use. Phenology analysis by using satellite time series data to assess regional scale vegetation variability and to identify change is consequently essential to characterize current patterns of land cover phenology and distinguish between anthropogenically driven land cover change, and interannual climatic fluctuations (de Beures & Henebry, 2004). Research in plant phenology has brought on intense interest in the context of climate change due to the relationship between the global climate change and vegetation phenology. The temporal and spatial characteristic of remote sensing data leads to this research topic because temporal and spatial characteristics are both important for the study of plant 11 phenology. Additionally, remote sensing technology has changed the observation of plant phenology from points (phenelogical observation station) to coverage (macro region) for better phenology observation across regions, countries, continents, and even across the globe (Zhang, 2006). The scientific community has been continuously interested in developing methods of employing satellite data to monitor and extract changes in LSP. Time series of vegetation indices (VI) based on satellite observations were first created in the 1970s and have continuously grown in importance for ecological and other biosphere-related research (Beck et el., 2007). In the 1990s, several types of methodological research presented phenological parameters derived NDVI (The Normalized Difference Vegetation Index) from AVHRR sensor. Myneni et al. (1997) used AVHRR-NDVI data to show a lengthening of the growing season relating to warmer air temperature over the northern latitudes from 1981 to 1991. Tucker et al. (2001) extended the study of the growing season trends to the period 1982-1999. He found that the vegetation index in higher northern latitude was responding to warmer temperatures causing the growing season to start earlier and continue longer. Reed et al. (1994) have derived quantitative phenological parameters from 1989 to 1992 over the United States. NDVI data from AVHRR (1 km) presented the onset of greenness, time of peak NDVI, maximum NDVI, rate of green up, rate of senescence and integrated NDVI. The results of this study had implications for large area land use mapping and monitoring. Moulin’s research (1997) was the first effort to assess the main phenological stages of the vegetation at the global scale. This study proposed a method to derive the start, the maximum, the end, and the length of the vegetation cycle by using AVHRR sensors. In the 2000s, satellite imagery available globally with frequently repeating cycles has provided data to examine and monitor phenological events over large regions (Reed et al., 2003). An 12 initial analysis of the MODIS-NDVI and EVI performance in the global scale was presented at this time. Between 2000 and 2003, Zhang et al. (2005) explored the response of vegetation phenology to precipitation across Africa using MODIS and rain fall data. The result demonstrated that vegetation phenology was strongly dependent on the seasonality of precipitation. This characteristic suggested that climate change may have a significant influence on vegetation phenology in this continent. Zhang et al. (2006) found that phenological parameters exhibited strong correspondence with temperature patterns in mid and high latitude climates, with rainfall seasonality in seasonally dry climates, and with cropping patterns in agricultural areas. Multi-temporal time series data have been used to study vegetation phenology in SEA. Sakamoto et al. (2006) estimated the spatial distribution of heading date and rice-cropping system employed in the Mekong Delta relative to seasonal changes in water resources in 2002 and 2003 using MODIS-EVI data. Xiao et al. (2006b) developed a new geospatial database of paddy rice agriculture for 13 countries on South and SEA using a time series of MODIS-derived vegetation indices, 500 m spatial resolution in 2002. Tottrup et al. (2007) presented a new approach to mapping fractional forest cover across the highlands of Mainland SEA using regression tree modeling and multi-temporal MODIS 250 m data. Although there are research related vegetation phenology in SEA, most of these studies focused on rice cropping system or forest areas. They have not addressed trends and change of vegetation phenology in regional scale as well as have not analyzed the relationship between phenology and climate variability, which is the important problems in this region. 13 2.1.2 Remote Sensing Approaches Although AVHRR data with frequent temporal coverage (12 hours, 1 week, 15 days, etc.) have been used extensively for monitoring vegetation; the coarse scale of these data (more than 1 km) with five spectral bands result in difficulties of assessing spatial variations in vegetation amount and condition (Huete et al., 2002). The AVHRR sensor was originally designed for weather and climate study, and aims to provide radiance data for investigation of clouds, snow and ice extent, temperature of radiating surface and sea surface temperature (Xiao et al., 2009). Therefore, AVHRR data have some limitations for vegetation studies such as lack of calibration, poor geometry, and high level of noise due to large pixel size and limited cloud screening (Xiao et al., 2009). In addition to AVHRR, the Vegetation Instrument (VGT) Sensors of SPOT4 and 5 and LANDSAT can provide phenological characteristics. However, the spatial resolution of VGT is coarse (1 km) and LANDSAT data (30 m) are limited for their temporal resolution, availability, and cost as well as the high spatial resolution (30 m) are too large in file size to be utilized in a regional scale and long term study (Wardlow et al., 2007). In contrast, MODIS provides 8-day and 16-day dataset at spatial resolution of 250 m with substantially geometric and radiometric properties relative to AVHRR data. Furthermore, MODIS data have been atmospherically corrected and screened for clouds. MODIS, launched in December 1999, has become the standard sensor for phenological studies. The MODIS sensor has 36 spectral bands, seven of which are designed for the study of vegetation and land surfaces (Didan & Huete, 2006). The MODIS NDVI and EVI products are computed from atmospherically corrected bi-directional surface reflectances that have been masked for water, clouds, heavy aerosols, and cloud shadows (LP DAAC, 2010). MODIS also includes a new EVI that minimizes 14 canopy background variations and maintains sensitivity over dense vegetation conditions as well as uses the blue band to remove residual atmosphere contamination caused by smoke and subpixel thin clouds (LP DAAC, 2010). MODIS data, therefore, provide an improved basis for monitoring global ecosystem dynamics (Zhang et al., 2006). Huete et al. (2002) compared the MODIS-NDVI and AVHRR-NDVI, they revealed that the MODIS VI products provided better results for biophysical information for land surface characterization. Not only were NDVI used for phenological monitoring, Xiao et al. (2006a) also indicated that the MODIS-EVI data provide a good result to detect leaf phenology of tropical forests in a moist tropical region in South America. For these reasons, MODIS-VI products are suitable to use in tropical region such as SEA. 2.1.3 Vegetation Index To monitor global vegetation conditions and display phenological changes, a vegetation index (VI) is necessary for this research. VI data can be used for characterizing land surface biophysical properties and processes, including primary production and land cover conversion. A VI is spectral transformations of two or more bands designed to enhance the contribution of vegetation properties (Huete et al., 2002). A VI is computed from combinations of visible red and near-infrared spectral measurements. The advantages of using these numerical transforms rather than the original spectral observations include the following: minimizing soil and other background effects, reducing data dimensionality, providing a degree of standardization for comparison, and enhancing the vegetation signal (Malingreau, 1989 as cited in Reed et al., 1994). 15 One of the more commonly used vegetation indices, NDVI, has been frequently used to evaluate phenological characteristics over large areas. Although NDVI is widely used for the study of vegetation and phenology dynamics, NDVI is not suitable for the tropical zone. This is due to the variation of atmospheric conditions with aerosols and clouds, and it is more sensitive to soil background (Kobayashi & dye, 2005 as cited in Xiao, 2006a). Xiao et al. (2009) also indicated the limitation and problem of NDVI for LSP, the saturation, in particular for tropical forests and in high biomass region with large values of leaf area index. EVI is an 'optimized' index designed to enhance the vegetation signal with improved sensitivity in high biomass regions and reduce atmospheric effects and soil background. EVI includes the blue reflectance band to correct for the influence of atmospheric aerosols on red reflectance (Huete et al., 2002). EVI is based on this equation (Huete et al., 2002): (2.1) where ρ is atmospherically corrected or partially atmosphere corrected (Rayleigh and ozone absorption) surface reflectance, L is the canopy background adjustment that addresses nonlinear, differential NIR and red radiant transfer through a canopy, and C1, C2 are the coefficients of the aerosol resistance term, which uses the blue band to correct for aerosol influences in the red band. The coefficients adopted in the EVI algorithm are, L=1, C1=6, C2=7.5, and G (gain factor) = 2.5. In addition, new MODIS products with collection 5.0 provide a better continuity with the standard EVI index, and a new index called EVI2 (Didan & Huete, 2006). The EVI index is replaced by the EVI2 index over cloudy pixel, snow/ice covered pixels in order to reduce the effects of the blue band over bright targets. This is because the blue band, as well as 16 the red and NIR bands exhibit proportionality problems over bright targets that are aggravated by the atmosphere correction resulting in abnormally high EVI values (Didan & Huete, 2006). This new product is based filtering scheme and a modified compositing method to deal with residual and mislabeled clouds. EVI2 is used this following equation (Didan & Huete, 2006): (2.2) The EVI profile, which has potential to represent the seasonal and annual dynamics, can be generated to present the total trends and changes of vegetation conditions. Xiao et al. (2006a) mentioned that the EVI data provided a good result in detecting phenology in a moist tropical region in South America by using time-series data of EVI from MODIS in 2002. Huete et al. (2002) compared the MODIS-NDVI and EVI temporal profiles. The results indicated that the EVI values exhibited a smoother, more symmetrical seasonal vegetation profile with a narrower, well-defined peak greenness period. In addition, EVI is less sensitive to residual atmospheric contamination due to aerosols from extensive fire (Xiao et al., 2009). For these reasons, MODIS EVI is the best choice for this research because it reduces soil and atmospheric effects as well as designs to improve sensitivity in high biomass regions and to reduce the canopy background signal (Huete et al., 2002). Additionally, EVI of MODIS has proved that it can help to reduce common problems with vegetation indices, particularly reducing smoke of biomass burning in tropical areas, while it is sensitive to canopy variations in high biomass regions. These characteristics are appropriate for studying a complex annual cycle with several growing seasons during the year in tropical zones. Therefore, MODISEVI data, particularly MOD13Q1 product (16-day composite data at 250 m spatial resolution) are appropriate to detect phenological characteristics in SEA. 17 2.1.4 Phenological Information Derived from Satellite Time-series Data 1) Growing Season Profile and Phenological Parameters To extract meaningful phenology information about the vegetation growing season, it is necessary to identify the appropriate phenological parameters. Myneni et al. (1997) and Tucker et al. (2001) were among the first to use satellite data to identify some of the key phenological parameters from NDVI time series. They proposed a set of parameters, such as the time of onset of green-up and senescence, maximum rates of green-up and senescence, and the amplitude of maximum NDVI. Zhang et al. (2003) proposed to use four key transition dates— green up, maturity, senescence, and dormancy—to monitor global vegetation phenology in the Northeastern United States. Although this method can provide vegetation dynamic over large areas, it focuses only on identifying phenology behavior characteristics. To more fully understand phenology patterns of each season requires more phenological parameters. Therefore, to extract seasonal parameters, it is important to understand the seasonal curve of phenology. Jönsson & Eklundh (2004 and 2010) clearly explained this growing season profile (Figure 2.1). The start of a season, marked by (1), is defined from the filtered or fitted functions as the point in time for which the value of the left edge has increased to a user defined level (a certain of the seasonal amplitude), generally set to 10% of the distance between the left minimum level and the maximum level. The end of the season (2) is defined by the time for which the right edge has decreased to a user defined level measured from the right minimum level. The length of the growing season (3) is the time from the start to the end of the season. The base level (4) is given as the average of the left and right minimum values. The middle of a season is obtained as the mean value of time, the position (5) between the position (8) and (9), for which the left edge has increased to the 80% level and the right edge has decreased to the 18 80% level, respectively. The maximum of fitted data (6) is the largest data value for the fitted function during the season; it may occur at a different time from (5). The amplitude (7) of the season is obtained as the difference between the maximum value and the base level. The left derivative (8) is the rate of increase at the beginning of the season; it is calculated as the ratio of the difference between the left 20% and 80% levels and the corresponding time difference. The right derivative (9) is the rate of decrease at the end of the season; it is calculated as the absolute value of the ratio of the difference between the right 20% and 80% levels and the corresponding time difference. The large integral (10) is the integral of the function describing the season from the season start to the season end; it is given by the area of the region between the fitted function and the zero level, represents the total vegetation production. The small integral (11) is given by the area of the region between the fitted function and the base level or the integral difference between the function describing the season and the base level from the season start to the season end. The small integral represents the seasonally active vegetation, which may be fairly small for evergreen areas. Therefore, in evergreen areas, the first integral may be small even if the total vegetation is large. 2) Data Smoothing and Filtering VI is normally composited over a set period of time to process growing season and to extract phenological parameters. These composited data often have errors from cloud and atmosphere. Due to these errors from cloud and atmospheric contaminations, and varying sun and sensor angles, time-series data are required to remove outlier and spikes and minimize the effects of anomalous values (Figure 2.2). 19 1) 2) 3) 4) 5) 6) Start of season End of season Length of season Base level Middle of season Maximum of fitted data 7) Amplitude 8) Left derivative 9) Right derivative 10) Large integral 11) Small integral Figure 2.1 A simple NDVI profile for a typical patch of vegetation Source: Jönsson & Eklundh, 2010 Figure 2.2 The results of data smoothing applied to EVI time-series, the start and end of the seasons are marked with filled circles; Source: Jönsson & Eklundh, 2010 20 A variety of methods have been developed to reduce satellite noise of time-series data, from simple linear smoothing window methods to more complicated analytical curve function methods. Fourier Transform (FT) has become one of the principal methods of phenological analysis because phenologies typically have strong seasonal cycles, particularly those comprised of monthly averaged samples from multiple years (Hermance et al., 2007). However, FT methods may be problematic, when applied to irregular or asymmetric VI data. This is because they depend critically on symmetric sine and cosine functions (Chen et al., 2004). In addition, Fourier analysis fails to characterize each annual VI trajectory separately; it can generate spurious oscillations in the VI time-series (Hird & McDermid, 2009). Moreover, the observed data in classical harmonic methods need to be sampled at uniform intervals of time, which often poses a problem for typical VI data, particularly when there are extensive data gaps (Hermance et al., 2007). A high-order spline fitting model applies the recursive least squares procedure for the initial inter-annual fit, and then uses the average annual curve as a baseline and asymmetrically weights all data points above and below the average in order to fit the upper envelope of the data (Bradley et al., 2007). Although exponential weighting in this approach is appropriate for fitting rapid green-up, one drawback is that it also fits data errors (spikes) that fall above the time series (Bradley et al., 2007). Moreover, this approach is considered computationally intensive and time consuming (Hermance et al., 2007). Asymmetric Gaussian Functions (AG) uses semilocal methods. Local model functions are fitted to data at intervals around maxima and minima in the time series and used to build a global function that describes the annual cycles of vegetation (Jönsson & Eklundh, 2010). The main advantage of this approach is that this function can be applied to data not uniformly 21 sampled in time, i.e., the spacing between time-stamps does not affect the output of the filter (Bachoo & Archibald, 2007). Beck et al. (2006) mentioned that AG is appropriate for describing vegetation dynamics at high latitudes. However, this method may be difficult to identify a reasonable and consistent set of maxima and minima to which the local functions can be fitted, especially for noisy data or for data from areas where there is no clear seasonality in land cover time series (Chen et al., 2004). The Double Logistic Function (DL) performs based on the minimum and maximum VI, two inflection points, and parameters related to the rate of increase or decrease in VI (Beck et al., 2007). DL, like AG, can deal with negatively-biased noise and work well to approximate winter conditions (Hird & McDermid, 2009). Therefore, this technique is appropriate for describing vegetation dynamics at high latitudes (Beck et al., 2006). The drawback of this technique are similar to AG in that it relies on the linear combination of local and independent intra-annual functions; in some cases, it fails to match the global waveform of numerous time series (Carrao et al., 2010). Another drawback is that this technique maintains the upper envelop of values, it is not suitable to apply on EVI, which have less negatively-biased noise and more erroneous spikes than NDVI (Hird & McDermid, 2009). Savitzky-Golay (SG) Filtering is one of the potential functions for generating fitted curve of time-series data. SG method performs least-squares polynomial regression on each point (local polynomial function) to determine the smoothed value and adapt for the upper envelop of the VI time-series profile (Jönsson & Eklundh, 2010). The SG filter applies an iterative weighted moving average filter to time series, with weighting given as a polynomial of a particular degree (Chen et al., 2004). 22 The main advantage of this approach is that it tends to preserve features of the distribution and capture rapid phenology changes, such as relative maximum, minimum, and width, which are usually flattened by other adjacent averaging techniques (Jönsson & Eklundh, 2004). More importantly, SG is provided to model the shape of VI profile, even in problematic environments with large variation in the magnitude and timing of vegetation growth. It can distinguish the signal from the seasonal variation and noise to reconstruct a clear and clean timeseries by linear combination of nearby values. When considering influence of phenological parameters and biogeographical regions for smoothing VI time-series data, SG has proved to be an essential method to deal with these factors. Not only does it preserve features of the distribution and capture rapid phenology changes, SG provided good results in the monsoon climate of SEA with multiple annual seasons during a year (Tottrup et al. 2007). In tropical zones, there are large numbers of low NDVI values due to atmospheric perturbations; SG effectively replaces these low values with higher NDVI values that are more typical of the cropland found in these areas (Chen et al., 2004). Bachoo & Archibald (2007) compared two different fitted models, SG and AG, for smoothing and extracting phenological parameters in southern Africa. The result exhibited that SG approach worked better than AG. Tottrup et al. (2007) claimed that SG method provided the good results in smoothing data and extract phenological parameters for mapping fractional forest cover in SEA. They also confirmed that the SG approach is applicable in this monsoon climate with multiple annual seasons during a year. However, this method is very sensitive to the size of window, that is, a wide window limits the ability of the filter to follow rapid but relevant changes in VI (e.g., a rapid 23 green-up), while a narrow window may cause the result to over-fit the time series and retain more noise (Hird & McDermid, 2009). Therefore, the size of window has to be applied with caution, and it should be tested before running the overall data (Chen et al., 2004). Based on SG mentioned above, this method is robust to reconstruct high-quality VI time-series data in this research. SG filtering can preserve features of the distribution and capture rapid phenology changes, and deal with large variation in the magnitude and timing of vegetation growth. Most importantly, it is proven to provide good results in the monsoon climate of SEA with multiple annual seasons during a year. Therefore, this method is selected to produce a smoothed curve, while capture rapid phenological change in SEA. 3) Extraction of Phenological Parameters Phenological parameters have been derived through various approaches from satellite time-series data after smoothing time-series profiles. The inflection point on a fitted curve is one of these techniques, which was developed by Badwar (1984 as cited in Reed et al., 2003). This technique defined the start of the growing season (SGS) as where the (left) time derivative transitions from 0 to a positive number. Similarly, the end of the growing season (EGS) is defined as where the (right) time derivative transitions from negative to zero (Reed et al., 2003; Moulin et al., 1997). Zhang et al. (2001, 2003) defined the phenological parameters: green up, maturity, senescence, and dormancy from the inflection point technique. However, the problem of this technique is the difficulty in extracting VI in evergreen, snow effects, and the slow rate of senescence in some biomes (Reed et al., 2003). Curve derivation phenology is also used to identify phenological parameters by examining curve characteristics where VI data exhibit a rapid sustained increase (Reed et al., 2003). White et al. (1997) and Moulin et al. (1997) defined SGS as the timing of the highest 24 positive derivative in the NDVI. In a similar way, the lowest negative derivative is considered the EGS. However, this method suffers a major drawback in that it is not known whether the rapid change is the result of the natural variability of the data or a significant change. In addition, when the VI signal fails to follow an abrupt and rapid increase or decrease, this method is difficult to determine SGS and EGS (de Beurs & Henebry, 2010). Threshold-base phenology is the effective technique to extract phenological parameters. This technique uses a pre-defined or a relative reference value for defining SGS and EGS (Reed et al., 2003). The SGS is defined as the day of the year (DOY) that the VI crosses the threshold in upward direction; whereas, the EGS is defined as the DOY that the VI crosses the same threshold in downward direction (de Beurs & Henebry, 2010). Various methods of threshold based technique have been developed to extract phenological parameters, for example, thresholds based on long-term mean VI, thresholds based on a baseline year, threshold based on reference value, and thresholds based on VI ratios (de Beurs & Henebry, 2010). This research applied thresholds based on VI ratios for studying temporal variations of phenological dates by finding the appropriate values for SEA. This technique computes ratio of the seasonal amplitude measured from the left and right minimum values to define the SGS and EGS respectively (Jönsson & Eklundh, 2010). The rapid growth is more important than the first leaf occurrence for this technique and this method can reduce the effect of soil background at lower vegetation signals. In addition, users can select appropriate threshold depending on spatial local climate and VI characteristics; therefore, this technique is flexible and can be applied in various regions (Jönsson & Eklundh, 2004). Moreover, thresholds based on VI ratios are considered as a consistent determination, independent of the geographic location and land cover of the observed study area (de Beurs & Henebry, 2010). 25 4) Trend Analysis for Satellite Time-series Data The trend analysis of the VI time series is important for monitoring subtle and long-term vegetation changes to understand the inter-annual changes of vegetation dynamics. Changing in vegetation phenology reflects land cover change that covers both drastic and slight changes on vegetation. The drastic changes refer to land cover conversions (deforestation, urbanization) and the slight changes are modifications of land cover without changing its overall classification (Martínez & Gilabert, 2009). A wide variety of procedures have been developed for trend analysis of timeseries data. The characteristics of VI time series data are defined as non-stationary, which present different frequency components, i.e., seasonal variations, long-term and short-term fluctuations, the patterns are typically seasonality, trends and abrupt changes or discontinuities resulting from disturbance events (de Beurs & Henebry, 2005). With these characteristics, simple linear regression is not appropriate for time-series data because it always results in parameter estimates but these parameters are not always significantly different from zero and the slope parameter fails to report the associated error and the overall error of the linear regression model (de Beurs & Henery, 2005). In addition, regression requires normality of the error distribution as well as constant variance of the errors and the method is very sensitive to extreme values (Qin, 2011). In contrast to linear regression, Mann-Kendall trend test (MK) is one of the efficiently non-parametric trend tests and it is more suitable for non-normally distributed data and censored data, which are frequently encountered in hydro-meteorological times series (Yue et al, 2002). Additionally, this statistic is considered to be the most useful trend analysis statistics for environmental time series (Luus, 2009). The rank-based nonparametric MK statistical test has since been widely applied to assess the significance of trends in hydrological and climatological 26 time series including water quality, stream flow, temperature, and precipitation (Latifovic & Pouliot, 2007; Xu et al., 2007; Yue et al., 2002; Yue et al., 2003). Recently, MK has been applied to VI time-series analysis (Julien & Sobrino, 2009; Martínez & Gilabert, 2009; Meng et al., 2011). The basic principle of the MK test is to examine the signs of all pairwise differences of the observed values, ranked in chronological order (Hirsch et al., 1982; Hirsch & Slack, 1984). This test determines whether decreases and increases in a dataset are significant according to the sum of sign differences between observation pairs (Luus, 2009). Since the test uses the relative magnitude of data (based on rank order) and not the actual data values, this approach can deal with outliers, missing values, non-normality, seasonality, and values below detection limits (de Beurs & Henery, 2005; Ngwenya, 2006). Therefore, this technique ignores the degree of the change, and takes into account only the yearly positive or negative change (Minnesota State University, 2009). Numerous studies have shown that the MK method is robust for trend and change analysis of VI image time series. Julien and Sobrino (2009) retrieved local land surface phenology trends for the whole globe. The MK trend tests in this study have been performed for time series of phenological parameters derived by NDVI, and trends were estimated by linear regression. Martínez and Gilabert (2009) indicated that the MK method is easy to calculate, robust against non-normality, and insensitive to missing values. The inter-annual series has been used to identify the trend of vegetation dynamics from NDVI time series, which is related to land-cover changes in Spain. This study applied the MK method for the strength and direction of a trend and the Sen’s method was used to determine the slope of changes. 27 In addition, MK was used to distinguish the characteristics of changes in climate variability and NDVI during the period from 1982 to 2000 in China (Meng et al., 2011). This study mentioned that the MK method is suitable to detect one abrupt change in a sample time series of anomalies that are independent and Gaussian distributed. Cui et al. (2012) applied the MK test to analyze vegetation change treads in different vegetation-climatic regions from 1982 to 2006 in Inner Mongolia by using the long time series NOAA/AVHRR NDVI dataset and vegetation-climate regions map. Trend analysis of NDVI and land surface temperature parameters from AVHRR was present in the application of the yearly land-cover dynamics methodology throughout the 1981-2001 periods (Julien et al., 2011). These time-series data were tested for trends using the MK trend tests. Non-parametric MK trend test was also applied to analyze monotonic greening and browning trends from global NDVI time-series from 1981 to 2006 (Jong et al., 2011). The results showed that NDVI time-series data were characterized by outliers, seasonality and serial auto-correlation; however, using the MK method took account of these effects. Due to various advantages of this technique, this research used MK to identify the trend of satellite time-series data (EVI, phenology, and rainfall seasonality). 2.2 Phenology and Climate Variability 2.2.1 Relationship between Phenology and Climate Vegetation phenology provides evidence not only ecosystem dynamics, biodiversity, and land use change but also play an important role in global climate change. Vegetation phenological events are effective and sensitive indicators for detecting and monitoring global environmental change driven mainly by temperature and precipitation (Zhang et al., 2005). A variety of studies demonstrated a clearly identifiable phenological response to global warming in high northern latitudes where temperature is normally the most important 28 climatic factor limiting plant photosynthesis (Myneni et al., 1997; Tucker et al., 2001; Robeson, 2004; Schwartz et al., 2006). These studies showed changes in temperature that are relevant for different phases of plant development. The timing of spring events (budburst, flowering) has been earlier since 1970s mainly due to global warming. For example, in Northern Hemisphere temperate land areas, the Spring Indices showed that first leaf dates and last frost dates were 1.2 and 1.5 days per decade earlier, respectively from 1955 to 2002 (Schwartz et al., 2006). This situation has accelerated spring and led to a longer growing season (Robeson, 2004). These characteristics are important indicators of ecosystem dynamics and global environmental changes. Although several research identified that the global temperature increase affects ecosystem in the higher northern latitudes (Myneni et al., 1997; Tucker et al., 2001; Yu et al., 2003), few studies have investigated the effects and implications of climate change on the lower latitudes where both temperature and precipitation play important roles in limiting plant biological processes (Yu et al., 2003; Zhang et al., 2005). Vegetation phenology in arid and semiarid as well as tropical climates is primarily controlled by water availability (Zhang et al., 2005). These ecosystems are dependent on rainfall to trigger the emergence of green leaves and control vegetation growth duration (Kramer et al., 2000). In addition, precipitation patterns in the rainy season are also critical to crop germination, growth, and harvest (Zhang et al., 2005). More importantly, extreme climate events (floods and droughts) and the extension of rainfall deficits or excesses also influence phenology in these regions (Zhang et al., 2005). According to Cleland et al. (2007), in tropical ecosystems, phenology might be less sensitive to temperature and photoperiod, and more tuned to seasonal shifts in precipitation. Such shifts are expected to occur in relation to rising global temperatures, but both the direction 29 and magnitude of change vary regionally. Corlett & Lafrankie (1998) mentioned that in most of tropical Asia, seasonality in rainfall is the major ecological factor. Furthermore, in the monsoon tropics, seasonal changes such as drought are severe enough to have a significant impact on plant growth. SEA is in the monsoon tropic region and so it is highly sensitive to climatic fluctuations, particularly rainfall changes causing floods and droughts. According to ADB (2009) and FAO (2007), most of SEA suffers from recurrent droughts and floods. Floods occur with greater frequency than ever before and damage agricultural areas in the river delta in this region, i.e., the Mekong, the Red, and the Chaopraya River Delta (ADB, 2009). In addition, droughts also affect agricultural areas in this region. The unavailability of real-time climatic and environmental data makes it difficult to monitor phenology response to climate change in SEA. Additionally, numerous studies have focused on the study of relationship between phenology and rainfall in high northern latitudes (e.g., North America, China) or in Africa; vegetation dynamics and rainfall seasonality in SEA have not been addressed. It is consequently necessary to monitor and assess the relationship between phenology and rainfall for a better understanding of environmental changes in this region. 2.2.2 Satellite Time-series Data for Vegetation Phenology Response to Rainfall Change Satellite time-series data enable researchers to investigate the response of vegetation phenology to climatic variation on a regional scale. NDVI, which is the most commonly used remote sensing derived measurement, has been linked to rainfall seasonality in numerous studies. Nicholson et al. (1990) demonstrated a strong statistical association between 30 rainfall and AVHRR-NDVI in West African Sahel and in East Africa on climatic time and space in 1982-1985. However, they found that NDVI exhibited a saturation response to rainfall— NDVI continued to increase as rainfall increased; after a certain threshold value, NDVI remained constant, while rainfall was increasing. Lee et al. (2002) and Yu et al. (2003) examined annual climatic variation, and its influence on plant phenology of the steppes of Inner Mongolia using AVHRR-NDVI. This research found that the major factors driving the changes in plant phenology is climate variation (precipitation and temperature). This was also demonstrated by Li & Xiao – Ling (2007), they used AVHRR-NDVI (from 1992 to 2001), precipitation, and temperature data to analyze the annual and monthly changes of vegetation in China. Both precipitation and temperature had different effects depending on different types of landform; moreover, precipitation dominated the changes of vegetation cover in high mountain areas. While several studies used rainfall data from meteorological stations to examine the response of phenology to rainfall, Zhang et al. (2005) demonstrated that MODIS-EVI and rainfall data retrieved from TRMM (Tropical Rainfall Measuring Mission) with 1x1 degree grid cells can be used to calculate both the onset and the end of rainy seasons. Phenological patterns were linked to the onset and the end of rainy seasons to assess the phenological response to precipitation at continental scales in Africa (Zhang et al., 2005). Phenological patterns from MODIS data and rainy season from TRMM were proven useful tools for monitoring the response of vegetation to precipitation. TRMM is a satellite jointly launched by U.S.A. and Japan in 1997 designed to observe and understand the tropical rainfall (NASA, 2011). For data quality, TRMM will reprocess all products with improved algorithms approximately once per year and use ground 31 based validation (NASA, 2011). This satellite dataset for rainfall has minimized the need of ground based measurement because in many parts of the world there are few rainfall observations or stations; the data are hard to obtain; and in some countries, war or economic crises lead to multi-year breaks in the records. My research used MODIS-EVI time series data and TRMM (TRMM-3B42-daily product) to study the response of phenology to rainfall variation in SEA. This product is TRMMadjusted merged-infrared (IR) precipitation and root-mean-square (RMS) precipitation-error estimates. The algorithm of this product uses TRMM merged high quality (HQ)/infrared (IR) precipitation with other satellites (e.g., Meteosat-7, NOAA-12). The daily product is accumulated from 3-hourly product (3B42-3-hourly) with a 0.25 by 0.25-degree spatial resolution (NASA, 2011). By coupling satellite-derived vegetation phenology and precipitation, the response of phenology to seasonal variation in rainfall is possible. Although several studies explored the response of vegetation phenology to precipitation, small temporal extent and coarse resolution of phenology and rainfall data limits the scope with regard to vegetation dynamics. Few studies examined vegetation dynamics and rainfall seasonality in tropical climate as in SEA. Rainfall is definitely a dominant independent variable to vegetation and the relationship with phenology can also help to determine the sensitivity of various vegetation formations to climate variability. Therefore, the study of ecosystems in SEA will be helpful in understanding the climate of SEA and how it affects the vegetation important in light of global climate change. 32 2.3 Environmental Model: The DNDC Model 2.3.1 Global Climate Change and Greenhouse Gases Emissions In the context of global climate change, increasing greenhouse gases (GHG) is the main cause of climate change (Li et al., 1994; Pathak et al., 2005). Several trace gases, such as carbon dioxide (CO2), methane (CH4), nitrous oxide (N2O), are the key greenhouse gases. CO2, CH4, and N2O contribute towards global warming at 60%, 15%, and 5% respectively, and concentrations of these gases in the atmosphere are increasing at 0.4%, 3.0%, and 0.22% per year, respectively (Pathak et al., 2005). Agro-ecosystems are considered to be the major source of the increase in GHG emissions into the atmosphere (Jagadeesh Babu et al., 2005). Agriculture contributes 92%, 65%, and 26 % of the total anthropogenic emissions of N20, CH4, and CO2, respectively (Zhang et al., 2002). CO2, CH4, and N2O are products of the biogeochemical cycles of carbon (C) and nitrogen (N) in agro-ecosystems and are produced in soils through decomposition, nitrification/denitrification, and methanogenesis (Li et al., 2004). In addition, any changes in management or climate/soil conditions will alter the biochemical or geochemical processes leading to changes in the gas fluxes (Li et al., 2004). Conversion of naturally vegetated land to agriculture generally results in the loss of soil organic carbon (SOC). The loss of carbon from agricultural soils contributes to the growth of atmospheric CO2 concentration and the potential for global warming (Li et al, 2003). In addition, food production contributes approximately 70% of global atmospheric input of N2O (Li et al., 2005). Variation in N2O fluxes is related to SOC content and C-sequestration during the 33 nitrification/denitrification process. Increasing C-sequestration is directly linked with increasing N2O emissions (Li et al., 2005). In particular, using fertilizer will alter the nitrification/denitrification process when SOC increases and results in increases of N2O (Li et al., 2005). Furthermore, several studies reported that agricultural activities are responsible for approximately 50% of global atmospheric inputs of CH4, and rice paddies are the major sources of CH4, contributing about 12% to global methane emissions (Zhang, et al., 2012; Zhang et al., 2009; Jagadeesh Babu et al., 2005). Additionally, numerous observations indicated that agricultural production and anthropogenic activities, e.g., tillage, fertilization, irrigation, and manure amendment, have resulted in significant increases in the number of GHG emissions (Li et al., 2004; Pathak et al., 2005; Zhang, et al., 2009; Zhang et al., 2012). CO2, CH4, and N2O are consequently the major GHG emissions from agricultural activities. In order to predict the rate of GHG emissions from various ecosystems, a number of environmental models have been developed in recent years to understand environmental consequences. The different modeling approaches can be grouped into the following categories: regression, empirical/semi-empirical, and process-based models (Cai et al, 2003). Early models used regression relationships between rates of emission and either crop biomass or grain yield. A simple regression model with a GIS framework was applied by Sozanska et al. (2002) to make an inventory of N2O emissions in British soils. The DAYCENT ecosystem model has also been applied to simulate soil organic carbon levels, crop yields, and annual trace gas fluxes for various soils (Del Grosso et al., 2002). CENTURY was developed for carbon, nitrogen, sulphur, and phosphorus cycles (Parton et al., 1988). Matthews et al. (2000) developed a process-based 34 Methane Emissions from Rice EcoSystems (MERES) model for simulating CH4 emissions from rice fields. Cao et al. (1995) estimated CH4 emission from rice fields in China and globally by using a process-based Methane Emission Model (MEM). Mall and Agarwal (2002) compared the performance of crop simulation models, CERES-Rice and ORYZA1N in India to predict the influence of climate change on rice grain yield. Recently, a process-based model INFOCROP has been developed for scaling-up gas emission estimates from tropical agriculture (Agarwal et al., 2004). 2.3.2 DNDC Model and its Capacities One of the most effective environmental models is the DenitrificationDecomposition Model (DNDC). This model is a process-based model developed by C. Li and his colleagues. The model is based on the biogeochemical concepts for predicting soil biogeochemistry (Li et al., 1992, 1994). The DNDC model was originally developed for predicting carbon sequestration, nitrogen dynamics, and trace gas emissions for agricultural lands, and then it was expanded to simulate N2O, CH4, CO2, nitrogen oxide (NO), and ammonia (NH3) emissions (Li, 2000; Li et al., 1992, 1994). This model is widely used to quantify GHG emissions from agriculture and other land uses. As a process-based biogeochemical model, the DNDC can simulate GHG dynamics under a range of changing ecological drivers (soil physical properties, climate, topography, vegetation) and management variables (cropping, tillage, manure, fertilization, and grazing practices), while capturing temporal and spatial variability (Li et al., 1992, 1994, 1996). The DNDC is able to estimate the rates of trace gas production and consumption in agricultural ecosystems by simulating the fundamental processes controlling the interactions among 35 ecological drivers, soil environmental factors, and relevant biochemical or geochemical reactions (Zhang et al., 2002). During the past two decades, the DNDC model was tested and validated by many researchers and under a wide range of conditions worldwide for predicting GHG emissions. The model has been applied to estimate soil organic carbon dynamics (Li et al., 1994, 2003), N2O emissions from agricultural fields (Li et al., 1996, 2001) and dairy farms (Brown et al., 2001), CH4 emissions from rice fields (Li et al., 2004,2005b; Pathak et al., 2005; Jagadeesh Babu et al., 2005; Cai et al., 2003; Fumoto et al., 2008; Smakgahn et al., 2009; Zhang et al., 2009; Zhang et al., 2012), and a forest version of DNDC, PnET-N-DNDC, was developed for simulating N2O and NO emissions from forest soils (Butterbach-Bahl et al., 2004). With continuous modification and calibration, the DNDC model can become powerful tool for applying in national trace gas inventory studies under different conditions and different areas, for example in the U.S., Canada, the United Kingdom, Germany, Italy, New Zealand, China, Japan, and Thailand (Li et al., 1994; Li et al., 2004, 2005a, 2005b; Cai et al., 2003; Butterbach-Bahl et al., 2004; Jagadeesh Babu et al., 2005; Pathak el al., 2005; Fumoto et al., 2008; Smakgahn et al., 2009). The DNDC model has been modified for predicting GHG emissions from paddy rice ecosystems (Li et al., 2004, 2005b). The new DNDC model modified simulations of anaerobic biogeochemistry and rice growth as well as the parameterization of paddy rice management. This modified model was tested for its sensitivities to management alternatives and variations in climate and soil properties, particularly in tropical areas. The modified DNDC model was developed and used for estimating emissions of CO2, CH4, and N2O from all rice 36 paddies in China, Japan, Thailand, and India (Li et al., 2004, 2005b; Cai et al., 2003; Pathak el al., 2005: Fumoto et al., 2007; Smakgahn et al., 2009; Zhang et al., 2009, 2011, 2012). The DNDC provides significant advantages over other environmental models for quantifying GHG emissions from agricultural ecosystems. DNDC has been linked to a crop model to simulate SOC dynamics and emissions of several trace gases from both upland and wetland agricultural ecosystems (Li et al., 1992, 1994, 1996, 2004). DNDC simulates soil redox potential and CH4 emissions from saturated soils and CH4 uptake; whereas other models, such as DAYCENT, process CH4 uptake in non-saturated soils (Olander & Malin, 2010). Furthermore, this model can manipulate the quantity of irrigation with different types of irrigation and can be applied to a wide variety of crop types (Salas et al., 2008). Moreover, the DNDC model can work in site mode or regional mode and can be linked to a GIS database of climate, soil, vegetation, and farming practices (Cai et al., 2003; Jagadeesh Babu et al., 2005). The previously described capabilities of the DNDC model are critical for quantifying emissions in tropical zones that have varieties of crop types in both upland and wetland ecosystems. In addition, the extensive validation and applications worldwide demonstrate that DNDC is capable of capturing the basic patterns and magnitudes of SOC, modeling C and N dynamics, and quantifying trace gas emissions across a broad range of climatic zones, soil types, and management regimes. 2.3.3 GHG Emissions in SEA and Thailand SEA is considered highly vulnerable to consequences of climate change (ADB, 2009). With 563.1 million people, the population is rising almost 2% annually, compared with the global average of 1.4% (ADB, 2009). Land use change is one of the most dramatic changes in this region – particularly change from forest to agricultural areas. Additionally, high 37 concentrations of human activities and economic growth in this region are highly dependent on agriculture, natural resources, and forestry. Accordingly, the agricultural sector is very important for providing jobs and livelihoods for a great number of people in this region. SEA contributes significantly to global GHG emissions due to the increase in large scale agriculture. According to ADB (2009) and IPCC (2007), SEA contributed 12% of the world’s GHG emissions in 2000, an increase of 27% from 1990, a rate greater than that of the global average. The land use change and forestry sector has been the major source of emissions from this region, contributing 75% of total regional GHG emissions in 2000 ADB (2009). Increasing demand for food and industrial crops has led to increases in agricultural land conversion and intensification, generating considerable environmental pressure. Nonagricultural areas (for example, forestlands, grasslands, and wetlands) have been converted to cropland for the production of beans, natural rubber, palm oil, rice (paddy), sesame seed, soybeans, and vegetables. This change results in emissions of CO2 from biomass (below ground and above ground) and soils. The burning of biomass in the process of land conversion also causes emissions of CO2 and other gases such as CH4, N2O, carbon monoxide (CO), and oxides of nitrogen (NOx). In addition, conversion of natural ecosystems to agriculture has caused changes in SOC and the decrease in SOC elevates increasing atmospheric CO2 (Li et al., 1994). SEA countries have to use modern crop varieties and improved farming techniques (e.g., irrigation, tillage, fertilizers) to increase yields. These processes have significant impacts on the environment, particularly CO2 and CH4 levels. 38 In addition to CO2 emission from agricultural areas, paddy rice fields are also an important source of CH4 (Li et al., 2005b). Rice is a major food crop in SEA. Most rice farms in Asia are flooded paddy fields, which are the major source of CH4 (Li et al., 2005b), and rice farm management, e.g., water management and fertilizer, has greatly affected the total CH 4 emissions (Li et al., 2004). Additionally, the greater use of nitrogen fertilizer or chemicals without efficient management will have significantly adverse impacts in N2O emissions (IPCC, 2007; ADB, 2009). Intensified cropping systems (i.e., two or three crops per year) with increased tillage, irrigation, and N fertilizer use have been widely applied in SEA to increase crop production. Such techniques can increase production, but they have contributed to decreases in SOC and to increases in the emissions of CH4 and N2O. These gaseous emissions not only represent potential economic losses, but also could lead to a negative impact on environmental quality (Li et al., 2004). This situation is particularly important to consider in Thailand because it is an agricultural country and rice production is its major agricultural product. According to Thailand Report (USDA, 2010; Towprayoon, 2005), paddy rice cropland in Thailand accounts for 52% of all cultivated land in the country (approximately 10 million hectares of rice-growing areas) and for 6% of the world’s rice paddies. Additionally, Thailand is the world's largest exporter of rice. The total rice production in 2010 was 30.7 million tons and the export was estimated at 8.8 million tons. Therefore, rice production in Thailand represents a significant portion of the Thai economy, which uses over half of its farmable land area and labor force, and rice is the main sources of nutrition for most Thai citizens. Most of rice growing areas are in the northeast of 39 Thailand, although the Central Plains is known as the nation's "rice bowl". Wetland (rainfed and irrigated rice fields) contributed about 99% of the total rice area. Irrigated rice fields are about 25% of rice lands, while the rest is still rainfed rice (FAO, 2007). These paddy rice croplands are considered to be the large CH4 source in Thailand. As one of the world's leading rice producers and one of 150 nations that have signed the United Nations Framework Convention on Climate Change in 1992, Thailand needs to estimate emissions at a broad scale in order to set mitigation policies and develop effective monitoring methods (UNFCC, 2007). Several approaches have been developed for estimating GHG emissions from cropping systems in Thailand, particularly CH4 from rice paddies. This is because wetland rice soils have been identified as an important CH4 source at the global scale (Smakgahn et al., 2009). Flooded rice fields, with their abundant organic matter, warm temperatures, and anaerobic conditions, provide and ideal environment for methanogenic activity (Matthews, et al., 2000). In addition, CH4 provides an equivalent warming effect being about 32 times higher than CO2 (Matthews, et al., 2000). Therefore, accurately estimating CH4 emissions from rice paddies has become important for GHG inventories or mitigation policies at country and regional levels. However, estimates of GHG emissions from cropping systems are still far from reliable because of large spatial and temporal variations of the emission records they are based on (Cai et al., 2003). Additionally, Smakgahn et al. (2009) also mentioned that CH4 estimations are difficult due to a variety of factors, including soil properties, rice varieties, and agricultural practices. Several studies attempted to conduct the estimation of GHG emissions from Thai rice 40 fields. A typical approach is to use emission factors and scaling factors for various categories of ecosystems provided by the guidelines of the IPCC to estimate CH4 emissions from rice fields in Thailand (Smakgahn, 1999; Gale et al., 2005; Eiumnoh et al., 2002). In addition, a variety of studies quantified CH4 from rice fields in Thailand based on experiments conducted in specific sites or other field methods to study the interaction of various factors and their effects on CH4 emissions in order to find the appropriate technology to reduce the emission of this gas (Jermsawatdipong, 2002; Tawprayoon, 2005; Siriratririya, 2001) Great efforts have been made to estimate CH4 emissions from Thai rice fields using the DNDC model at Prachinburi, Suphan Buri, Chiang Mai, and Surin (Cai et al., 2003). Smakgahn (2003) also simulated CH4 emissions using DNDC model at Kanchanaburi, Nonthaburi, Samutsakorn, and Singburi. However, the results showed large discrepancies between the simulated results and filed observations. The simulation of emissions using the DNDC model was revised and improved in the research of Smakgahn et al. (2009) to estimate CH4 from irrigated rice fields at nine study sites in Thailand and the results showed a high correlation with observations. The results suggested the considerable uncertainties resulted from model description (root biomass) and input parameters from soil properties and field management (reducible soil iron, straw incorporation rate, and water regime). In addition, scientists have measured CH4 and N2O from Thai rice fields by using process bases and empirical models (Towprayoon & Smakgahn, 2003; Smakgahn, 2003). The recent research estimated CH4 emissions by using vegetation indices from LANDSAT-5TM (Keereerom, 2011). 41 Spatio-temporal analysis of emissions from rice fields with changing environmental conditions and field management at regional and national scales has become important for GHG inventory and mitigation and for understanding environmental change in Thailand. However, most research has focused on quantifying emissions at a local scale (Jermsawatdipong, 2002; Cai et al., 2003; Tawprayoon, 2005; Smakgahn et al., 2009). The studies that used a model for GHG estimation frequently mentioned the problems of model accuracy due to missing data, incomplete databases, and insufficient availability of input parameters for the model (Smakgahn, 2003; Smakgahn et al., 2009). The local circumstances such as soil types, climate data, and farming management are necessary to better simulate emissions but these data are unavailable in some sites or some areas. In addition, a large spatial dimension is required to meet the demands for GHG mitigation. Remotely sensed data have provided an effective way to derive regional input data for the DNDC model. Spatial and temporal information obtained from satellite data provide geographic characteristics such as rice growing areas, rice phenology, and soil moisture, particularly in places where field observations are unavailable. Many researchers have used remotely sensed data for mapping the extent of paddy rice on local and regional scales and applied these data to the DNDC model (Zhang et al., 2009; Zhang et al, 2011; Zhang et al., 2012). In addition to rice areas, phenology is very important data for DNDC to quantify emissions and these data are difficult to access or are often unavailable. With spatial and multitemporal remote sensing data, phenological patterns (e.g., the start and end of the growing season) can be extracted at regional scales and this information provides the opportunity for the DNDC model to generate spatio-temporal patterns of emissions at regional scales. 42 Therefore, this research applied remote sensing technology to develop the new database based on a grid-based system (250 m x 250 m) for the DNDC model in order to identify spatial patterns and changes in GHG emissions. This approach is different from other methods in terms of scale. Most DNDC modeling performed at site and regional scales usually is applied to a county as the basic unit. Spatial variation in climate, soil, and farming management provides great uncertainties at this scale unit, for example, using average soil properties for each county. In this way, the simulation at the site level with a grid-based unit can provide more detailed and accurate agroecosystem information because this approach could reflect the spatial diversity of crop growth environments. This approach is able to reserve the advantages of site-based modeling but also meet the demand for large-region estimation. In conclusion, agroecosystems play an important role in balancing food production and environmental protection. As a result, predicting and estimating the impacts of trace gas emission are significant for environmental safety and decision making. An environmental model is required to estimate regional GHG emissions in agricultural areas, particularly in Thailand. The new approaches proposed in this research, with spatial and multitemporal remote sensing data and the new database system with a grid-based unit at 250 m resolution can improve model performance to estimate inter-annual variation in emissions from rice paddy fields in Thailand or elsewhere. Additionally, the distributions and changes of these GHG emissions can be investigated with phenological changes and climate data in order to study the effects of the environmental impacts and ecosystem dynamics in this region. 43 2.4 Visualization for Satellite Phenological Information 2.4.1 Remote Sensing Time-series Data and Visualization Satellite time-series data are important for regional, natural resource management and environmental monitoring because they can represent changes in space, time, and attribute. These data show processes, patterns, trends, directions of change, anomalies, and causes and relationship of dynamic phenomena on the earth’s surface (Harrower, 2002; Blok, 2006). To date, there are no more technological or data problems. The important problem is how to represent time-series data to show the distribution and changes of geographic phenomena. To represent the spatio-temporal processes of environmental systems, cartographers need to use dynamic geovisual displays to realistically mimic the passage of world time. This goevisual display is important for environmental monitoring and management because it provides the capability to visualize time-series data that congruently depict change over time and space (Harrower, 2002; Blok, 2006). For example, map animation allows cartographers to create displays of time-series data that mimic this dynamic process, whereas static maps of timeseries data have the limitation to show change over time. Additionally, dynamic geovisual displays enable viewers not only to detect changes but also to understand the meanings encoded in transitions within these dynamic displays (Fish et al., 2011). Although satellite time-series data can dynamically display spatio-temporal processes in the real world, there are difficulties for the viewer to understand the meanings encoded within dynamic displays. The major challenges of satellite time-series data for visualization are data characteristics that require users’ knowledge and experience (Blok, 2006). In addition, satellite time-series data also require effective techniques to design the visualization (Harrower, 2002). Phenological information derived from satellite time-series data is the useful 44 information, but it does not has easily available method for users to view, explore, and share the data. To respond this need, it is necessary to develop the effective methods to visualize these geographical phenomena of the changes. 2.4.2 Visual Exploration Visual exploration is an effective technique for time-series satellite images. Blok (2006) described that visual exploration is a creative process to derive meaning and construct knowledge. Visual exploration is a process that users start without much knowledge about the underlying data, but interactivity provides the tools for the search of structures and trends. The comparison of patterns, spatial and temporal characteristics, relationships and trends, enables users to interact with the data and their graphic representations. Visual exploration is potentially helpful in all stages of the knowledge construction process (DiBiase et al., 1992). Users can derive the meaning of information from this process to formulate hypotheses and can also be used to confirm, synthesize, and ultimately present ideas and information in the research process. Due to now the rate of environmental change effected by human activity has been increasing, we need more data to depict and monitor the changes. However, quantitative data about this change exceed our capacity to learn from them. Using visual exploration to convert these data to graphic form and searching visually for patterns and anomalies can increase potential for scientific insight in the earth science (DiBiase et al, 1992). Harrower’s research (2002) demonstrated that the effects of change on multitemporal remote sensing data enabled users to detect structures or trends and helped users to formulate specific hypotheses about the behavior of geographic entities represented within the data. Therefore, satellite time-series data have potential for visual exploration to detect the patterns of changes. 45 2.4.3 Cartographic Design for Satellite Time-series Data In order to create visual exploration, the geographical representation in Internet atlas with map animation and bivariate mapping make it possible to detect spatially or temporally distributed regional patterns, and changes of satellite time-series data. Map animation and bivariate mapping are powerful techniques applied to the interactive atlas because they can display dynamically spatio-temporal processes and allow viewers to easily understand complex spatio-temporal geographic phenomena. This research creates an Internet atlas combining thematic maps, bivariate maps, and map animation and distributing this information by web mapping application with the MapServer. Five componets (Internet atlas, thematic maps, bivariate maping, map animation, and web mapping application) are discussed in the following sections. 1) Internet Atlas The development of Internet atlas is one of the cartographic developments on the web. An atlas is a combination of maps and additional information combined with a well structured work (Richard, 1999). The electronic form of the atlas applying the techniques of multimedia is commonly referred to as a multimedia atlas or Internet atlas (Borchert, 1999). An online atlas is a collection of maps, illustrations, informative tables or textual matters bundled with a clear structure and connected on an HTML page (Richard, 1999). The main advantages of the Internet atlas are data and media exploration by selection of context (spatial or thematic), interactive functions and dynamic animation for dynamic contents or method comparisons, individual selection of geometry and attributes, and individual map design by users (Borchert, 1999). Due to these advantages, the production and uses of an Internet atlas are currently more available on the Internet for the general public. 46 At present, there are several atlases that can be found on the web and they are presented in an attractive way with useful information. The classification of Internet atlas/multimedia atlas can be distinguished by view-only atlases and interactive atlases. The interactive atlas is used for analysis purposes, which sometimes include a mapping module, GIS functionality, and topographic base map collections. This research created thematic map, bivariate map, and map animation in the Internet atlas with interactive functions for users to explore phenological information and related environmental variables in SEA. 2) Thematic Map The thematic map is used to display the spatial pattern of a theme or attribute and can provide information in three basis ways: specific information, general information, and pattern comparisons on two or more maps (Slocum et al., 2009). The thematic map was applied in this research to present one specific data set, emphasizing the spatial pattern of one or more geographic attributes such as EVI, start and end date of the growing season. These maps are appropriate to indicate the pattern and magnitude differences of geographic variable. According to Slocum et al. (2009), there are the particular concerns in thematic mapping. The first concern is the classification techniques. Thematic mapping is needed to standardize data and select appropriate techniques for data classification. The appropriate technique depends on data characteristics and purpose of map. The other concerns are an appropriate number of classes, level of measurement, visual variables, including the color scheme. For example, hue is appropriate for qualitative attributes but lightness and satuation suggest quantitative differences. The last issue is legends that should be designed so that users can easily interpret maps. These numerous factors should be carefully considered, when designing the map. 47 3) Bivariate Map Slocum et al. (2009) defined bivariate mapping as the process of displaying two attributes of geographic phenomena. The purpose of bivariate mapping is to compare and examine the relation between two attributes in order to represent individual distributions or the correlation between them. Consequently, bivariate maps enable the visualization of statistical correlation and show a lot of information on one map. Bivariate mapping is the technique that was developed by the U.S. Bureau of the Census in the 1970s (Slocum et al., 2009). The graphical representation of bivariate data can enhance the ability of map readers to detect and comprehend important phenomena as well as communicate major conclusions to others (Wang, 1998). However, bivariate mapping has been criticized because of its failure to communicate information, although it can be effectively displayed with a clear legend and explanatory note (Olson, 1981). The disadvantages of this type of map are the subjective interpretation, distortion of the data representation, and erroneous impression (Olson, 1981; Wang, 1998). In addition, it is difficult to read and to design, particularly when the number of entities to be represented increase; moreover, a lot of the correlation depends on classification methods (Slocum et al., 2009). To deal with the significant challenges of bivariate mapping, cartographers should concern with cartographic design for this type of map. Olson (1981) proposed the early revised design for bivariate maps: a very prominent and clear legend, the displays of single-variable and two-variable maps on the same layout to increase and improve extraction of information (singlevariable maps should show in black and white colors), and explanatory notes to reinforce the confidence of the user and show the types of information extractable from the display. 48 The important issues for bivariate mapping are data classification and color schemes (Brewer, 1994; Slocum et al., 2009). An appropriate method of classification should be selected for bivariate mapping. Furthermore, no more than three classes are appropriate for bivariate maps because it can increase cognitive overload leading the users to be confused with the maps. Color schemes for bivariate maps should suggest correlation and the technique is dependent on types of map and types of data. Brewer (1994) introduced color schemes for bivariate choropleth maps focusing specifically on correlation between variables. According to Brewer’s research (1994), a comprehensive set of color scheme types and corresponding guidelines for the use of hue and lightness for each scheme are proposed such as sequential/sequential schemes and diverging/diverging schemes. Brewer (1994) concluded that different characteristics of distributions and their interrelationships can be examined by using different schemes. In addition to data classification and color schemes, the important issues of bivariate mapping are the selection of appropriate variables and the use of a legend (Slocum et al., 2009). Both issues are dependent on the specific purpose and map-users. With these effective designs, bivariate map can communicate the relation between variables. Furthermore, users can analyze the data, detect and comprehend important phenomena, and find the major conclusion for their research. 4) Map Animation To represent spatio-temporal processes, cartographers use map animation to communicate changes in space, time, and attribute simultaneously (Goldsberry & Battersby, 2009). Animations are defined as sequences of static graphic depictions (frames) to show the graphic content in rapid succession (Harrow, 2008). Animation provides a means for better 49 understanding the complexity of geographic changes because it can represent the outcomes of change and the process of change (what happened, how, and why), while the static map only represents the outcomes of change (Harrow, 2002). In addition, animated maps seem especially suited to emphasizing change over time and space, whereas static maps cannot represent time and geographic behaviors or processes. Although small multiple maps can represent changes over time, they can show this change only in macro steps (Harrower, 2008). The representation of small multiple maps is dependent on how many key events are selected to discretely display geographic phenomena. This characteristic is different from animated maps, which show continuous changes or micro steps in complex systems that might be missed in small multiple displays (Slocum et al., 2009). Furthermore, small multiple maps are difficult to interpret when comparing subregions within each map (Slocum et al., 2009). Animated maps are effective in conveying concepts that are dynamic, particularly when using satellite time-series data because represent geographical processes that change over time and space that are important to investigate the complex and dynamic behavior of geographic entities (Harrower, 2002). This dynamic process is important for many applications such as deforestation, erosion, and land use/land cover change (Blok, 2006). Additionally, map animation represents changes in phenomena of interest i.e., space, time, and attribute (Goldsberry & Battersby, 2009), for example, selecting year (time) to explore space and attribute changes. With these capabilities, map animation is significant for knowledge construction and the research process. It can be used to formulate specific hypotheses and confirm information in the research process. 50 Map animation can communicate effectively the visual dynamic of a geographic process; however, it is difficult to perceive and easy overload users with rapid changes. Several researchers have argued that effective map design can reduce these challenges (Harrow, 2002, 2008; Blok, 2006; Goldsberry & Battersby, 2009). For example, using the number of classes appropriately, applying dynamic visual variables and interactive map animations, and determining the levels of change detection can facilitate visualization in map animation. Additionally, understanding the meaning of map animation depends on the symbology used, the level of interaction that users are permitted to have with the animation, and the expertise and experience of the users in data characteristics and animation. These techniques can apply to represent changes of satellite time-series data with map animation and allows users to detect meaningful information. 5) Web Mapping Application: A Tool to Disseminate Geographic Data Web mapping is the process of creation, distribution, and use of web maps on the Internet. In addition, visualization aspects are the main subject of this mapping (Hachler, 2003). Web mapping opens the opportunity to disseminate Internet atlas, thematic map, map animation, and bivariate map to users with easy accessibility but cheap distribution of data and software. One type of web mapping application is a web map server. A map server makes available geographical information in the form of maps. These maps can be static maps as a result of a cartographic process or interactive maps that users can change and define some selected graphical variables of the map presentation (Gartner, 1999). The “MapServer” of the University of Minnesota (UMN) is a well-known open source project whose purpose is to display dynamic spatial maps over the Internet The UMN MapServer has been developed in cooperation with NASA and the Minnesota Department of Natural Resource to make its satellite 51 imagery available to the public (The University of Minnesota, 2011). This open source product has already been successfully used in several other projects including in environmental monitoring and management tasks (Hachler, 2003; Choi et al., 2005; Aden et al., 2010). The main advantages are to operate on UNIX/LINUX and Windows systems and provide core functionality to support a variety of web mapping applications. The most important features include vector and raster format support, fully customizable and template driven output, automatic legend, and scale bar building (Hachler, 2003). While the growth rate of web mapping and map servers have been increasing, there is no web mapping of phenological information, especially in SEA. With powerful spatial data representation and simple user interfaces, web mapping can increase the understanding of spatial phenomena and facilitate environmental monitoring and management (Tsou et al., 2004). Therefore, the web mapping application with the MapServer that facilitate exploration of satellite time-series data and present valuable phenological information will benefit potential users to explore and analyze the temporal and spatial change of geographic phenomena. In addition, users can specify their unique areas and variables of interest and then see the information described in a map with easy-to-navigate tools and simple user interfaces. In conclusion, satellite time-series data that represent vegetation conditions and phenological parameters are suitable to create the specific visualization system. This visualization system provides opportunity to explore geographic phenomena, particularly the distribution and pattern of changes. Visual exploration enables users to derive the meaningful information from interactive processes by the comparison of patterns, and their spatial and temporal characteristics, as well as relationships and trends of the data. Internet atlas is one of the most effective systems, which provides data exploration for analysis purposes with applicable 52 functions. In addition, this system is opened to general public by using the MapServer to deliver spatial information on Internet. Additionally, thematic mapping, map animation, and bivariate mapping can enhance the capability of internet atlas by enabling users to understand the complex spatio-temporal geographic phenomena more easily. With this effective system, users not only detect the trends and changes of environment (e.g., vegetation condition and phenological parameters) from satellite time-series data but also understand the meanings within these dynamic displays. Better understanding in environmental change by using visualization is the critical need to solve environmental problems and to provide effective planning and decisionmaking (Harrower, 2008). 53 CHAPTER 3 RESEARCH METHODOLOGY There are four main parts in this research: spatio-temporal variation of vegetation phenology in SEA, driving forces of phenological changes, environmental impact, and visualization system. The flow chart of research methodology is shown in Figure 3.1. The first part of this research explored spatio-temporal variation of EVI and phenological parameters; then, causes and drivers of phenological changes were examined by addressing the relationship between phenology and rainfall seasonality and investigating phenology and land use changes. Next, environmental model, DNDC, was used to assess environmental impact in rice fields in Thailand. The change and distribution of GHG emissions were examined with phenology and climate variation. The last part is the visualization system. Interactive Phenological Atlas for SEA (IPA) was established to display the outputs from previous sections to present spatiotemporal patterns of environmental changes in SEA. 3.1 Study Area Southeast Asia (SEA) consists of two dissimilar portions: the Indochina Peninsula and the Insular SEA (Archipelagic Nations). This research focuses on the seasonal dynamic in the Indochina Peninsula of SEA: Thailand, Vietnam, Cambodia, Lao PDR, Myanmar, and Malaysia; (Figure 3.2). This is because phenology can be extracted in the Indochina Peninsula while there is no variation in seasonal vegetation cycles in the Insular SEA due to the rainy tropical climate along with rain forests. For convenience, those countries in Indochina Peninsula are referred to as “Southeast Asia (SEA)" in this research. 54 Figure 3.1 Research methodology 55 Figure 3.2 Study Area, the Indochina Peninsula in SEA The Peninsula has a humid subtropical climate with a winter dry season and much of it receives considerable annual precipitation (Southeast Asia, 2009). With these characteristics, the Peninsula has tropical maritime climate featuring relatively high temperatures, high humidity, and abundant precipitation. There are three basic seasons in this region: the rainy season, winter, and summer. Most of SEA is covered with tropical forests. There are rain forests (tropical evergreen forests), which have a high annual rainfall, and monsoon forests (tropical deciduous forests) where rainfall is seasonal (Southeast Asia, 2009). 56 According to NIC (2009), SEA has a natural climate variability due to the regular pattern of seasonal monsoons or the periodic shift in global climate caused by ENSO (El Niño-Southern Oscillation events, a global climate phenomenon that recurs irregularly every 2-7 years and is associated with changes in sea surface temperature and prevailing winds). Seasonal monsoons can cause extreme weather events, such as floods and droughts. Moreover, ENSO can intensify existing floods and droughts in this region. 3.2 Data Sources The following data are used in this research (Table 3.1): Table 3.1 Summary of data requirements in this research 1 2 3 4 5 Data Description EVI images retrieved from MODIS dataset onboard NASA’s Terra spacecraft (MOD13Q1 product) (http://reverb.echo.nasa.gov/reverb/) Rainfall data retrieved from TRMM (TRMM3B42 daily product) (http://mirador.gsfc.nasa.gov/cgibin/mirador/presentNavigation.pl?project=TRMM &tree=project) Rainfall data acquired from meteorological stations (Thai Meteorological Department ) LANDSAT images (http://glovis.usgs.gov) Land use data in Thailand - National Land Use Dataset of Thailand (Land Development Department, Ministry of Agriculture) 57 Spatial and temporal coverage Indochina Peninsula: 250 m spatial resolution from 2001-2010, every 16day EVI composite period (23 dates per year with 460 tiles per year) Indochina Peninsula: daily rainfall rate for 0.25x0.25 degree spatial resolution from 2001-2010 Thailand: Monthly rainfall: 20012010 Thailand: 2004,2005,2006, and 2007 Thailand: 2001, 2002, 2009, and 2010 Table 3.1 (cont’d) 6 Data Description Data for DNDC Model - Soil property data (Department of Land Development) - Rainfall data (Thai Meteorological Department ) - Land use data (Land Development Department, Ministry of Agriculture) - Farming practices (e.g., tillage, fertilization, residues, flooding period) (Department of Agriculture, Department of Agricultural Extension, Rice Department, Local rice research center, Office of Agricultural Economics (Ministry of Agriculture and Cooperatives), research of Jermsawatdipong et al. (2002)) Spatial and temporal coverage Lopburi Province in Thailand Daily rainfall: 2001-2010 2002, 2010 2002, 2010 3.3 Data Pre-processing The flowchart of this step is shown in Figure 3.3: data pre-processing for EVI and phenological analysis, and data pre-processing for rainfall analysis. EVI time-series data of 4,600 tiles were downloaded and organized by dates and years, and collected in a database preparing for pre-processing. The first step of time-series EVI data processing was reprojection from the sinusoidal projection to the geographic projection, WGS84. The next step was creating mosaics from each tile to produce one image per 16 days (23 images per a year). The mosaic process produced 230 time-series files which were subsetted to the study area. After pre-processing, TIMESAT program was used to remove noise, generate EVI profiles, and extract phenological parameters. 58 Figure 3.3 Flowchart of phenology and rainfall analysis 59 3.4 Fitted Function and Phenological Extraction This research employed the TIMESAT program (Jönsson & Eklundh, 2010) to remove noise, generate EVI profiles, and extract phenological parameters. TIMESAT is an open source software developed by Jonsson and Eklundh, and was designed to process time-series of vegetation index derived from satellite spectral measurements (Jönsson & Eklundh, 2010). This program is a simple yet robust method to remove noise in time-series data. The major advantages of this program are a user-defined parameter in the smoothing process and extracting phenological parameters, providing three different smoothing functions to fit time-series data (asymmetric Gaussian, double logistic, and adaptive Savitzky–Golay filtering), and providing a comprehensive set of phenological parameters (Tan et al., 2011). In this research, the 16-day MODIS EVI composites were used to depict phenology and seasonal variations in vegetation activity. TIMESAT needs time-series data of three complete years to extract phenology; therefore, three years of EVI time series data were used to process, and a number of annual seasonal parameters were extracted from the middle year. For example, to extract phenology in 2005, we need EVI time-series data from 2004 to 2006. This research implemented separate three years of EVI time series to extract phenology in the middle year as shown in Figure 3.4 (e.g., using data from 2001 to 2003 to extract phenology in 2002 and using data from 2002 to 2004 to extract phenology in 2003). TIMESAT first implements the simple median filtering to remove noise, spikes, and irregular values of original images. Then, a filtering technique was used to produce a smoothed curve of EVI profiles. To obtain accurate EVI profiles and phenological parameters, the appropriate threshold parameters were required to be defined in the TIMESAT Program according to the local conditions in SEA before applying the smoothing filter to data series. A 60 number of settings were generated based on the local conditions in SEA by exploring time-series data and tuning and testing the important parameters. Several tests were run on specific areas in Thailand according to available data for assessing model performance. The parameters in the setting files were modified for SEA (Table 3.2), for example, window size of SG technique, the seasonality parameter for more than one growing season per year, amplitude cutoff value to remove pixels with very weak seasonality from the processing. 0.7 2001 0.6 2002 2003 EVI 0.5 0.4 0.3 0.2 0.1 0 1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 MODIS Date Figure 3.4 Three year time-series data to extract phenological parameters from one year 3.4.1 Data Filtering After obtaining appropriate setting for parameters, smoothing functions were processed for time-series data. TIMESAT was developed to make the fit of the upper envelope approach of the time-series and to reflect significant anomalies through an iterative process. This approach is more flexible and effective in obtaining high-quality time-series. The TIMESAT program uses three processing steps based on least-squares fits to the upper envelope of VI time series data as follow: - Removing Spikes and Outliers: As some spikes and outliers may be detected in time-series; these positive and negative outliers seriously impair the function fits. The 61 TIMESAT program implements the median filtering to remove these spikes and outliers. Data values that differ from the median value by more than the spike value multiplied with the standard deviation of VI values and that are different from the left and right neighbors are removed (Jönsson & Eklundh, 2010). Table 3.2 TIMESAT input parameters Parameters Job_name Image /series mode (1/0) Trend (1/0) Use mask data (1/0) Data file list/name Mask file list/name Image file type Byte order (1/0) File dimension (nrow ncol) Processing window (start row stop row start col stop col) No. years and no. points per year Valid data range (lower upper) Mask range 1 and weight Amplitude cutoff value Print functions and weights (1/0) Output files (1/0 1/0 1/0) Use land cover (1/0) Spike method Spike value No. of landcover classes Land cover code for class 1 %Season parameter No. of envelope iterations (1-3) Adaptation strength (1-10) Force minimum (1/0) and value Fitting method (1-3) Weight update method Window size for Sav-Gol. Season start method Season start / stop values 62 Values SEA 1 0 0 filelist2001.txt none 2 0 13473 9766 3 23 1 10000 0 1500 0 110 0 1 2 1 1 0.1 3 2 0-99999 1 1 4 1 0.2 0.2 - Determining the Number of Seasons: The number of growing season is determined by using de-trended data values (ti, yi), i = 1, 2, . . . , N for all years in the timeseries to fit to a model function by using the following equation (Jönsson & Eklundh, 2010): f(t) = c1 + c2 sin(ωt) + c3 cos(ωt) + c4 sin(2ωt) + c4 cos(2ωt) (3.1) Where, ω = 6¶/N, t = a function of time, de-trended data values = (ti, yi), i = 1, 2, . . . ,N for all years in the time-series. The first basis function determines the base level whereas the pairs of sine and cosine functions correspond respectively to one and two annual vegetation season equation (Jönsson & Eklundh, 2010). The fitting procedure always gives a primary maximum. In addition, a secondary maximum may be found, if the amplitude ratio between the secondary maximum and the primary maximum exceeds a user defined threshold; therefore, we have two annual seasons. - Savitzky-Golay Filtering: Fitted curve function of the Savitzky-Golay (SG) filtering is selected to remove noise, spikes and irregular value of original images due to cloud contamination, atmospheric conditions and bidirectional effects. The SG function in the TIMESAT program is suitable for studying vegetation dynamics. In contrast to functions resulting from Fourier methods, the SG function in the TIMESAT is able to capture interannual change, i.e., changes in seasonal timing in different years. In addition, the SG computes leastsquares polynomial regression on each point to determine the smoothed value and adapts for the upper envelope of the EVI time-series profile. According to the SG function in the TIMESAT program (Jönsson & Eklundh, 2010), the general mathematical process is: 63 ∑ (3.2) Where, Yi (new EVI values) is replaced by the average of EVI values in the window at position i, yi is the original EVI value at i position, i=1,…,N by a linear combination of nearby values in a window, N = the width of the moving window, determine the degree of smoothing, which consists of 2n+1 points, the weights are Cj = 1/(2n + 1). The moving average method preserves the area and mean position of a seasonal peak, but it does not maintain the width and height. In order to preserve these properties, SG performs least-squares polynomial regression on each point to determine the smoothed value for each point. For each point value yi, i = 1, 2, . . . , N we fit a quadratic polynomial f(t) = c1 + c2t 2 + c3t to all 2n + 1 points in the moving window and replace the value yi with the value of the polynomial at position ti. 3.4.2 Phenological Parameters Considering the profile of the growing season (Figure 2.1) and the review of phenological applications, this research identified 7 phenological parameters (Table 3.3) which are appropriate to depict phenological vegetation variability in SEA. This research applied a threshold-based function by calculating the ratio for the start (SGS) and the end of growing season (EGS) where the fitted curve reaches a proportion of the seasonal amplitude measured respectively from the left and right minimum value. The advantage of TIMESAT program is that it is easy to tune the threshold according to the local conditions. 64 Table 3.3 Phenological parameters Phenological Parameter 1. Start of the growing season Definition The first date on which the value has increased to a defined level measured from the left minimum level 2. End of the growing season The last date on which the value has decreased to a defined level measured from the right minimum level 3. Length of the growing season 4. Mid of the growing season 5. Amplitude of the growing season Time from the start to the end of the season The middle date between the start and the end date The difference between the peak value and the start and the end date 6. Large integral 7. Small integral The total vegetation production The seasonally active vegetation, the different between the season and the base level In addition, the most important parameter for extracting phenology is the range of the growing season. TIMESAT needs a specific date between which the season is expected to occur, and pixels that have the season in a given range will be processed. To ensure that the full season in SEA is captured, the range of the start and end dates should be wider than the expected season. There are multiple cropping frequencies in this region and the length of the second growing seasons is longer than one calendar year. For example, some vegetation may start in October and end in March the following year. The range of season should be large enough to allow for a certain variation in the start and end of the first and the second seasons over the processed area. Therefore, time series data were processed to test and tune as well as to compare with the land use data and phenological data in Thailand. After testing the range of growing season in SEA, the results show that a large number of pixels show no phenological data when the narrower range was applied (Figure 3.5). The range from January to July of the following years is the appropriate range of growing season to detect phenology and is large enough to cover the start and end dates of the growing season in this region. Due to the limitations of the 65 TIMESAT program, this research focused on the first and second seasons in this region but the third season cannot be extracted by TIMESAT program. Input settings for the TIMESAT program in SEA is shown in Table 3.2. Testing Phenology Range, Start Date, S1, 2006 Start date S1 Jan-Dec Jan-Dec Jan Jan-Dec Feb Mar April May June July Aug Sep Jan-Dec Jan-Dec Jan-Dec Oct Nov Dec Jan Feb Mar April Figure 3.5 Testing the range of growing season in Thailand 3.5 Rainfall Seasonality Extraction The rainfall profiles (daily rainfall from TRMM), which has potentials to represent the seasonal and annual dynamics, were generated to present the patterns and changes of rainfall seasonality from 2001 to 2010. First, this research extracted the mean annual rainfall to indicate the amount of rainfall throughout the season. Then, the start (SRS), the end (ERS), and the length (LRS) of rainy season were extracted. 66 Rainfall is a discrete event; therefore, the appropriate thresholds must be defined to derive seasonality. The criteria to define thresholds are that rainfall threshold should directly influence the vegetation seasonality (Zhang et al., 2005). Moreover, it should be an independent estimation of rainfall and should be dealt with sporadic events of rainfall, which can be misleading. To define thresholds for the rainfall seasonality in SEA, the criteria of biological method were applied in this research. This research applied the criteria of rainy season provided by Jutakorn (2011) to determine rainy season in SEA. Jutakorn investigated SRS, ERS, and LRS including variation of rainy season from 1951 to 2009 in Thailand. In this study, SRS is defined as the duration of which rainfall persists for five consecutive days and the amount of cumulative rainfall is higher than 10 mm, which is sufficient for vegetation to start the growing season. The first day of that rainfall period is the start of rainy season. This criteria identifies SRS from April, 1st of every year as the earliest month that rainfall could start in this region. ERS is defined under the similar criteria except that it is considered backwards, starting from December, November, etc. The last day of the duration of which rainfall persists for five consecutive days but the amount of cumulative rainfall is greater than 10 mm is the ERS. The period between SRS and ERS defines the length of rainy season. 3.6 Spatial Variation and Trend Analysis 3.6.1 Trend Analysis This research applied the Mann–Kendall (MK) trend test and the Sen’s method for EVI, phenology, and rainfall analysis. According to Yue et al. (2002), Xu et al. (2007), EPA (2006), and De Beurs et al. (2009), the MK trend test is calculated by listing the observations in temporal order, computing all differences that may be formed between measurements and earlier 67 measurements, and then counting the number of positive and negative values. The only information about these differences used in the MK calculations is the sign (positive or negative). If a value of the observation is higher than a previous observation, one is added to the test statistic and a plus sign is scored. If the values are equal, no sign is scored, and if the value is lower than a previous observation, one is subtracted from the test statistic and a minus sign is scored. Next, this method counts the number of positive and negative values. The sum of these values indicates the trend of the data. The positive sum indicates an increasing trend and the negative sum indicates the decreasing trend. MK is generated in the following form (Yue et al., 2002; de Beurs & Henery, 2005; Xu et al., 2007). The MK test is based on the test statistic S defined as: ∑ ∑ (3.3) where xj are the sequential data values, n is the length of the data set (number of years), and sgn (xj - xi) = +1 xj - xi >0 sgn (xj - xi) = 0 xj - xi =0 sgn (xj - xi) = -1 xj - xi <0 (3.4) The statistic S is approximately normally distributed with the mean and the variance as following: 68 [ ∑ ] (3.5) th T is the number of tiled groups, and tj is the number of points in the j group. These data are used to test the null hypothesis, Ho: there is no trend, HA: there are downward or upward trends at an α = 0.05 significance level. If |Z| > Z α /2, then reject the null hypothesis of no trend. The test statistic Z of a particular pixel is calculated by (EPA, 2006): Z= - (S) (3.6) √ ar(S) Where sign(S) = 1 if S >0, 0 if S = 0, and -1 if S < 0 As MK technique can provide the strength and direction of a trend but not the magnitude, the overall slope was obtained by using the Sen’s method to quantify the magnitude of the change, which can be positive or negative. The pixels with significant trend were used to calculate the magnitude of slope. According to the statistical analysis guidance for environmental data (EPA, 2006), Sen’s method is a nonparametric alternative for estimating a slope. This method calculates slopes for all the pairs of time points and then uses the median of these slopes as an estimate of the overall slope. As such, it is insensitive to outliers and can handle a moderate number of values below the detection limit and missing values. The assumptions of this method are that there are n time points, and let Xi denote th the data value for the i time point. If there are no missing data, there will be n(n-1)/2 possible pairs of time points (i, j) in which i < j. The slope for such a pair is (EPA, 2006): bij = (Xj – Xi) / (j - i) 69 (3.7) where bij = slope between data point Xi and Xj, Xi = data measurement at time i, Xj = data measurement at time j. Sen’s slope estimator is the median of the n(n-1)/2 pairwise slopes. If there is no underlying trend, there would be an approximately equal number of positive and negative slopes, and thus the median would be near zero. 3.6.2 Spatio-temporal Patterns, Changes and Trends of EVI Profile The EVI profile, which has potential to represent the seasonal and annual dynamics, was generated to present the total trends and changes of vegetation conditions in SEA from 2001 to 2010. The 23 fitted images per year were calculated to obtain the mean annual EVI by using TIMESAT program from the previous step so as to indicate the overall greenness of vegetation throughout the season and to illustrate the interannual variability for 10 years. The EVI changes and trends were identified by using the MK trend analysis and Sen’s method. The results indicate spatial distributions and interannual change of EVI for ten years. 3.6.3 Spatio-temporal Patterns, Changes and Trends of Phenological Parameters Seven phenological parameters were derived from the EVI time-series of 10 years by using TIMESAT. The mean values of each year in each parameter were computed to show spatial variation for 10 years and changes were processed by MK trend analysis and Sen’s method. The results exhibit spatial distributions and the trends of phenological changes. 3.6.4 Spatio-temporal Patterns, Changes and Trends of Rainfall Seasonality Spatial distribution of mean annual rainfall, SRS, ERS, and LRS, were generated for ten years. The rainfall seasonality was processed to detect trends by using MK and Sen’s method. Spatial distributions and the trend of rainfall seasonality were presented from this step. 70 3.6.5 Relationship between Phenology and Rainfall Seasonality This step used the coefficient of correlation (R) to calculate the (3.8) relationship between each pair of parameters (mean annual EVI-mean annual rainfall, SGS-SRS, EGS-ERS, and LGS-LRS). The coefficient of correlation uses this equation (EPA, 2006): ̅ ∑ (3.8) ̅ √∑ Where, ̅ and ̅ ∑ ̅ ̅ are the sample means of Xi and Yi In addition, the significance of relationship p-value (probability associated with significance) was performed to test the significance of the correlation at a level of 95%. Therefore, the pixels with a p-value less than 0.05 will show significant correlation. The tstatistic for the slope was calculated by the following equation (EPA, 2006): √ (3.9) 2 Where, r = Coefficient of correlation, r = Coefficient of determination, and n = Sample size Furthermore, linear regression was performed on each pair of parameters for pixels with significant correlation. The slope from linear regression identified the trend of relation in order to estimate rate of change of phenological parameters with respect to corresponding change in rainfall seasonality for a period of ten years. 3.6.6 Analyze the Drivers of Phenology Changes Two related factors (climate and land use) were investigated to analyze the drivers of changes in selected sites. Land use changes were compared with rainfall variability and phenological trends. Rainy seasons from TRMM were analyzed with phenolgical patterns. SRS, 71 ERS, LRS, and annual rainfall were processed to obtain rainfall patterns and then compare with phenolgical parameters. The correlation coefficients between these two variables were calculated. Land use changes were also considered in order to identify the relationship with phenological changes and to clearly explain the driver of changes, which are human management or climate. 3.6.7 Comparison between MODIS Phenology and Field Observations and TRMM Rainfall and Station Rainfall Although validation is a serious impediment of land surface phenology derived from satellite images over large areas due to scale and dates of phenological events (Zhang et al., 2003; Reed et al., 2009), many studies have employed various techniques to validate satellitederived phenology such as using climate data (de Beurs and Henebry, 2005) or ground observations (Beck et al., 2007; Fisher and Mustard, 2007; Kang et al., 2003). This research attempted to compare the phenological parameters derived from MODIS EVI with field data in order to observe errors and quantify the agreement between satellite derived phenology and field observations. For the comparison between MODIS phenology and field data, the field data were obtained by interviewing famers and head of villagers in central and northeast Thailand. The sites were selected based on land cover diversity and the phenological changes. Two hundred points of phenological data (SGS and EGS) were compiled from various crop types (e.g., rice, cassava, sugarcane, corn) in 2011 to identify the phenological data in 2010 (Figure 3.6). The 2 coefficient of determination (R ) was used to test the agreement between MODIS phenology and field data in 2010. This coefficient of determination quantified how much of the variance in the observed dates is explained by the model. 72 In addition to comparison between MODIS phenology and field data, this research also compares TRMM rainfall with station rainfall. Annual rainfall data from 2001 to 2010 were processed along with the average ten-year mean of rainfall seasonality (SRS, ERS, and LRS) for these comparisons. Monthly rainfall was obtained from Thai Meteorological Stations (123 stations cover the whole Thailand) and rainfall seasonality was acquired from the research of Thailand Meteorological Department (TMD) (Jutakorn, 2011). The comparison was assessed based on the performance of the coefficient of determination, showing the agreement between TRMM rainfall and station rainfall. Figure 3.6 Field points in central and northeastern Thailand 73 3.7 Modeling DNDC 3.7.1 Model Description This research applied DNDC model to estimate CO2, CH4, and N2O emission from rice fields in Thailand (Figure 3.7 and 3.8) because of its capabilities to quantify emission in tropical zones that have varieties of crop types in both upland and wetland ecosystems. DNDC consists of two major components (ecological drivers and soil environmental variables and soil environmental factors to trace gases) with six sub-models—soil climate, plant growth, decomposition, nitrification, denitrification, and fermentation (Li, 2000; Li et al., 1994, 2004; Zhang et al., 2006). The interaction of the six sub-models enables DNDC to simulate a relatively complete suite of biochemical and geochemical processes, which occurs under both aerobic and anaerobic conditions. 1. The first component in this model is to link between ecological drivers and soil environmental variables. DNDC integrates the ecological drivers in the three submodels: The soil climate submodel calculates soil temperature and moisture profiles based on soil physical properties, daily weather, and plant water use. The plant growth submodel tracks crop growth and partitioning of the biomass into grain, stalk, and roots. The decomposition submodel simulates decomposition of soil organic matter driven by the soil microbial respiration. The three submodels interact with one another to finally determine soil temperature, moisture, pH, redox potential (Eh), and substrate concentrations in the soil profile (e.g., ammonium, nitrate, dissolved organic carbon) at a daily time step based on ecological drivers (e.g., climate, soil, vegetation, and anthropogenic activity). 2. The second component of this model is to link soil environmental factors to trace gases with the three submodels. The nitrification submodel calculates growth of nitrifiers 74 and oxidation of ammonium to nitrate. The denitrification submodel operates at an hourly time step to simulate denitrification and the production of nitric oxide, nitrous oxide, and dinitrogen. The fermentation submodel simulates methane production and oxidation under anaerobic conditions. These three submodels predict nitric oxide, nitrous oxide, methane, and ammonia fluxes based on the environmental variables in the soil. DNDC allows users to define the farming management and other parameters for a certain crop. As a result, this model requires the following input parameters: soil properties (e.g., soil bulk density, pH, texture, and SOC), climate data (daily temperature and precipitation), farming practices (e.g., crop type and rotation, tillage, fertilization, manure amendment, planting and harvest dates, irrigation, flooding, grazing, and weeding). Model outputs are crop productivity (grain, stem and root yield, N-uptake, N-fixation), trace gas fluxes (NO, N2O, CH4, CO2, and NH3), soil organic C and N pools, soil inorganic N content (Li, 2000). DNDC model was modified to efficiently simulate C and N biogeochemical cycling in paddy rice systems by adding a series of anaerobic processes. Additionally, rice growth sub-model was developed to quantify important dynamic parameters for modeling gas production and oxidation (Zhang, et al., 2002; Li et al., 2004; Cai et al., 2003). In addition, this model modification improves the simulation for crop growth and soil processes and enhances its ability to estimate CH4 from rice paddy field under a wide range of climatic and agronomic conditions (Fumoto et al., 2008) 75 Figure 3.7 DNDC model 76 Figure 3.8 DNDC flow chart 77 3.7.2 Study Area Lopburi province was selected to quantify emissions. Lopburi is located in the Chaopraya River delta in the central plains of Thailand. This study site has been rice fields for years and these rice fields produce one of the highest quality rice in the country and even the world due to its location and farming practices. In addition, the representative site exhibits the changes of phenological and climate variations over a ten year study. Therefore, the patterns and trends of phenology and climate in this site indicate the changes of the growing season, which are useful for the interpretation in this research. Most of rice fields in this site are irrigated rice and double rice cropping systems, particularly in Amphoe Ban Mi and Tha Wung (Figure 3.9). Single rice cropping systems are found in the rest of the province, particularly in Amphoe Khok Samrong and Sa Bot. Modern cultivation managements are extensively practiced in this location. 3.7.3 Database Development and Input Parameters Due to several spatially explicit environmental and cultural variables, a process base model, like DNDC, needs a great number of data to estimate emissions in large scales. This is a major challenge for using ecosystem model at a regional scale. Furthermore, the input data for DNDC at a regional scale are large spatial and temporal variations with different factors, such as, soil properties, rice varieties, and agricultural practices (Cai et al., 2003; Smakgahn et al., 2009). These datasets are required to initialize and run the model. 78 Rice Cropping frequency in Lopburi 2002 Single Double Rice Rice Figure 3.9a Rice cropping frequency in Lopburi province 2002 79 Rice Cropping frequency in Lopburi 2010 Single Double Rice Rice Figure 3.9b Rice cropping frequency in Lopburi province 2010 80 This research integrated biogeochemical models with remote sensing technology to advance the DNDC regional application. DNDC model was used to identify spatio-temporal patterns of GHG emissions (CO2, CH4, and N2O) from rice fields in Thailand over the past 10 years (2001-2010). New method and database systems were developed in this research to increase efficient estimation of DNDC model. Satellite images are one of the effective methods to provide important data such as rice growing areas, rice phenology, and soil properties. In addition, satellite time series data could provide spatial and temporal information to extrapolate the understandings of spatial distribution in a large scale dimension. Spatial and temporal characteristics of phenology derived from remote sensing data were used in DNDC to quantify emissions. Also, this research took a new step to use regional database based on a grid-based system instead of using average input parameters in a large spatial unit (e.g., county unit). The different characteristics in a large spatial unit could be significant for biogeochemical modeling. Thus, the simulation of grid based system can differentiate the difference of soil properties, climate data, and farming management and substantially improve the accuracy of the GHG estimations from DNDC model at regional scale. Before simulating the whole dataset, this research tested and run model in specific locations of the study area to access model performance and test sensitivity of input parameters. In addition, this research applied model to estimate emissions for both single and double rice cropping systems. It is required to search for the technique running DNDC for two cropping systems, which the growing season is not completed in a year. These steps were performed in order to obtain accurate results from the model. - Construction of Database for Regional Simulation: This research used the DNDC model version 9.4 with the site mode (a site-scale simulation) to simulate emissions of 81 CO2, N2O, and CH4 from rice fields in Thailand. Figure 3.8 shows flow charts of database manipulations used in estimating GHG emission from rice fields in this research. Although regional mode in DNDC model can process large area estimation, each spatial unit in this mode is assigned to the same farming practices and this mode uses default farming practices instead of user specific crop management parameters (Xu, 2011). Consequently, this leads to the difficulty in adjusting models to specific sites. Therefore, this research simulated GHG emissions at the site level with grid-based unit to provide the output at regional level. This approach takes advantage of site-based modeling and also provides output for large-region estimation. DNDC used databases with spatially and temporally differentiated information on climate (temperature and precipitation), soil properties, and farming practices (crop type and rotation, tillage, fertilization, irrigation, and planting and harvest dates) as parameters for supporting regional scale analyses. The database was constructed based on a grid-based system. Each input data (e.g., phenology, soil, climate, farming management) was built based on grid unit (250 m x 250 m) as described in the following (Table 3.4, Figure 3.8). - Input Parameters: input data are described in the following: 1. The climate dataset composed of daily maximum and minimum air temperature and precipitation from 2001-2010 acquired from weather stations in Thailand (Table 3.1). The station weather data were interpolated to climate zone by using Thiessen (Voronoi) polygon analysis. Thiessen polygons are used to divide space into a number of regions according to a set of points. Thiessen proximal polygons are constructed by creating a triangulated irregular network (TIN) and generating the perpendicular bisectors for each triangle edge and forming the edges of the Thiessen polygons (ArcGIS Resource Center, 2011). As a result, Thiessen polygons contain only a single point input feature and any location within a Thiessen polygon is closer to 82 its associated point than any other point input features. The thiessen polygon method was used to transform point-based climate station to climate zone in this research. Each grid cell was assigned to use climate data from a weather station in its climate zone. In other words, grid cells within one polygon (one climate zone) will receive the same climate data. The interpolation was performed by using geoprocessing functions provided by ArcGIS 9.3. 2. Spatial soil databases, which consisted of soil profile, soil series, soil names, and soil physical and chemical properties (e.g., soil pH, texture, organic matter content), were induced from the 1:50,000 Thailand soil database. 3. The rice fields were obtained from National Land Use Dataset of Thailand (Table 3.1). Phenological information (planting and harvest dates) were derived from MODIS EVI time-series data (see more details in section 3.4.2) to map the spatial pattern of start and end dates of the growing season. 4. Agricultural management practices on rice cultivation and crop parameters of rice were investigated from Thailand reports and papers as well as by communicating with researchers, local agronomists, and local rice research centers in Thailand (Table 3.1). The detailed information on rice cultivation for DNDC model is shown in Table 3.4. In addition, annual yields were applied in crop parameter setting. However, unavailable input data were replaced by default values provided in the model. All the spatially differentiated input information was constructed in a grid-based system with a cell size of 250 m x 250 m. All datasets were segmented to rice areas and ignored other croplands. Then, the input text files were generated for DNDC model by using ERDAS and Python programing to derive data in each grid cell into text files. Therefore, one text file is for one site or one grid cell. For example, there are 40,000 cells in Lopburi province, which in turns, 83 generate 40,000 input text files. All input text files, including climate data, soil properties, and crop management parameters, were linked to DNDC through its interface with batch running for site mode (Figure 3.8). Table 3.4 Input parameters for DNDC model Parameters Environmental Factors Climate Soil Agricultural Management Rice area Rice cultivation (Planting/Harvest Dates) Tillage Fertilization (kg N/ha) Residues Management (Fraction of residues left in fields) Flooding a d Data Source 2002-2010: Daily climate data (Max-Min Temperature, rainfall) Rainfall data acquired from meteorological stations (Thai Meteorological Department ) pH, Density, SOC, % Clay (Department of Land Development) Land use data 2002, 2010 National Land Use Dataset of Thailand (Land Development Department, Ministry of Agriculture) Phenological Data 2002-2010 The results from the previous section a 15 DBP / plow depth of 20 cm 2002 b c #1 20 DAP ; AP = 90 #2 60 DAP; Urea = 30 2010 #1 20 DAP; AP = 150 #2 60 DAP; Urea = 60 10% - Local rice research center - Rice Department, Office of Agricultural Economics, Ministry of Agriculture and Cooperatives - Jermsawatdipong et al., 2002 Start: 7 DAP d End: 15 DBH Continuously flooding b c DBP = days before planting; DAP = days after planting; AP = Ammonium phosphate; DBH = days before harvest 84 The model simulated the emission cell-by-cell across the entire rice fields in the study site. All input files provided ecological, soil, and environmental drivers to process the crop and soil biochemistry sub-models of the DNDC model. The model, then, simulated C and N circulation of paddy rice system at a site scale for each grid cell and simulated the crop growth and environmental impacts of rice fields. The spatio-temporal patterns of CO2, N2O, and CH4 were processed using ARCGIS and Python programming to show maps of the annual emission rate per pixel (kg C/ha/year). The yearly total emissions (kg C/year) in each cell were calculated by multiplying the emissions by the size of a cell containing rice fields (multiply by the size of the grid cell if that cell has rice fields). The regional emissions of province were derived from a sum of all cells in the province for a simulated year. In addition, this research assessed Global Warming Potential (GWP) to indicate the net effect on global warming. GWP is defined as the cumulative radiative forcing integrated over a period of time from the emission of a unit mass of gas relative to some reference gas (IPCC 1996). In other words, GWP is the sum of the warming forces of all the three GHGs (CO2, N2O, and CH4) and convert them to a common basis of CO2-equivalents (Li et al., 2004). This method is used to estimate the warming effects of different long-lived greenhouse gases relative to each other within a particular time frame (for example, 100 years) as compared with the radiation a unit emission of CO2 absorbs in the atmosphere over the same time period. CH4 and N2O can be converted into CO2-equivalents according to their radiative forcings. The warming forces of CH4 and N2O are 25 and 298 times higher than that of CO2 per unit of weight with the 100-year time horizon (IPCC, 2007). The GWP value is 85 calculated by summing up the CO2-equivalents of all the three GHGs as shown in the following equation (adjusted from Li et al., 2004 according to the new standard of IPCC (2007)): GWPi = CO2 i /12 x 44 + N2O i/28 x 44 x 298 + CH4i /12 x 16 x 25, (3.10) where GWPi (kg CO2 equivalent/ha/yr) is the Global Warming Potential; CO2i, N2Oi, and CH4i are CO2 flux (kg C/ha/yr), N2O flux (kg N/ha/yr), and CH4 flux (kg CH4-C/ha/yr), respectively. 3.7.4 Spatio-temporal Patterns and Changes of GHG Emissions under Different Scenarios In order to identify the effect of the input parameters on the model regional emissions, this research tested five scenarios for a selected site (Table 3.5). The first scenario is the baseline scenario for 2002 and 2010, which was constructed based on the actual climate, phenology, land use, and farming management (fertilization, flooding period, and yield). The second scenario tested the effect of fertilization on GHG emissions by using different set of fertilization for 2002 and 2010 and keeping all other input parameters constant (apply input parameters of 2002 for 2010, except for fertilizer rate). The third scenario conducted multiyear simulations from 2002 to 2010 to investigate the impacts of interannual variations in climate condition on GHG emissions. This scenario used actual nine-year climate data (temperature and rainfall) for each simulated years, while keeping all other input parameters constant (used input data from 2002 for all years from 2002 to 2010). The fourth scenario was also set for multiyear simulations from 2002 to 2010. This scenario examined the effects of phenological changes by using actual phenology (start and end of the growing season) for nine years and all other inputs were held constant by using input data from 2002 for all years from 2002 to 2010. The last scenario investigated two parameters, phenology and climate for simulation of nine years. The 86 actual climate and phenological data for the nine-year period were utilized in this scenario and the remaining input parameters used were from the 2002 dataset. The results of these scenarios could represent the change of GHG emissions with the effect of human management as well as climate variation. Table 3.5 Five simulated GHG emissions from rice field in Thailand under different scenarios Parameters Scenario1 Baseline (2002, 2010) Climate Rice area Phenology (Planting/Harvest Dates) Fertilization Other management practices e.g., residues, flooding, tillage Conditions Scenario2 Fertilizer Change (2002, 2010) Scenario Scenario3 Climate Change (2002-2010) Scenario4 Phenology Change (2002-2010) 2002, 2010 2002, 2010 2002, 2010 2002 2002 2002 2002-2010 2002 2002 2002 2002 2002-2010 Scenario5 Phenology and Climate Change (2002-2010) 2002-2010 2002 2002-2010 2002, 2010 2002, 2010 2002, 2010 2002 2002 2002 2002 2002 2002 2002 Use actual data for 2002 and 2010 Use different rate of fertilizer (kg N/ha) 2002 #1 AP = 90 #2 Urea = 30 2010 #1 AP = 150 #2 Urea = 60 Multi-year simulation From 20022010 and apply actual climate data for each simulated years Multi-year simulation From 20022010 and apply actual phenology data for each simulated years Multi-year simulation From 20022010 and apply both actual climate and phenology data for each simulated years AP = Ammonium phosphate 87 3.7.5 Comparison of DNDC Results to IPCC Approach and Thailand Research Due to the unavailable field data for validation, this research assessed the results of DNDC model by using IPCC guidelines and the standard emission rates from the papers in Thailand. The modeled CH4 emissions were compared with IPCC method based on baseline emission factor (IPCC, 2006). The IPCC emission factor contains more specific values for CH4 emission in different types of rice cultivation. The result of this research is also compared with standard estimation factors for rice ecosystem in Thailand in 2005 and estimated CH4 emission rates for each province in Thailand based on minimum, median, and maximum scenarios in 1998 (Gale et al., 2005). It should be noted that few studies have been conducted on emissions of CO2 and N2O from Thai rice fields and no reliable estimations of their contribution to global emission are available. 3.8 The Interactive Phenological Atlas for SEA (IPA) 3.8.1 Conceptual Framework One objective of this research is to develop the prototype of Interactive Phenological Atlas for SEA—with facilitated, interactive, and usable interface to represent spatial distribution and changes in phenology and related environmental variables in SEA (Figure 3.10). This prototype is presented as Internet atlas supported by MapServer. The web prototype includes basic functions—zoom, pan, and simple selection—and advanced functions—map animation, dynamic legend, and graph creation. The graphic user interface of the web map provides functionality of legends in interactive web maps combining the controlling function for thematic map, bivariate map, and map animation. 88 Figure 3.10 Interactive Phenological Atlas for SEA The purpose of this visualization system is for visual exploration. Visual exploration is a creative process to derive meaning and construct knowledge (Blok, 2006). Users can derive the meaningful information from interactive processes from the comparison of patterns, spatial and temporal characteristics, relationships, and trends. This process enables users to formulate hypotheses and construct knowledge as well as to confirm information in the research process (DiBiase et al., 1992). In addition, this prototype is designed for the general public. The intended audience is environmental scientists and remote sensing and GIS specialists as well as students who would like to use these data for their works and research. 89 There are three data groups in this system (Table 3.6). The first group consists of thematic maps of EVI, phenology, rainfall, GHG emissions, trend maps, bivariate maps, and data quality assessment maps, which were implemented by web map server. The second group consists of map animation generated by using flash program. The third group is vector data (supporting data), Indochina boundary and Lopburi province. Table 3.6 Data and map contents in visualization system Data Details I. Maps 1. EVI images from 2001-2010 - Mean EVI of each year (2001-2010) and ten-year mean EVI 2. Phenological parameters - Spatial pattern/distribution of vegetation phenology from 2001-2010 (Start, end, length, mid, amplitude, large, and small integral) 3. Rainfall - Annual rainfall, the start, end, and length of rainfall seasonality (2001-2010) 4. GHG emissions - Spatial pattern and distribution of CO2, CH4, and N20 (Lopburi province, Thailand) - EVI trend (2001-2010) - Phenological trend (2001-2010) - Rainfall trend (2001-2010) 5. Trend maps 6. Bivariate maps - Start date and end date (2001-2010) - Length and amplitude (2001-2010) - Length and large integral (2001-2010) 7. Data Quality Assessment Maps - Data quality assessment (2001-2010) and ten-year mean data quality assessment II. Map animation - EVI and phenology (start date, end date, mid date, amplitude, length, large integral, small Integral) (20012010) III. Supporting data 5.Indochina boundary 6.Lopburi province 90 3.8.2 Cartographic Design for Satellite Time-series Data In order to serve maps in this visualization system, cartographic design is needed for thematic map, bivariate map, and map animation. - Thematic Map Design: Satellite time-series images in raster format were generated to create thematic maps for users to display on web map application. Thematic mapping can display the pattern and change over space. However, this mapping method needs appropriate techniques to design. The factors involved in thematic mapping design are data classification techniques, number of classes, color scheme, and legend design. The visualization system in this research generated thematic maps of EVI, phenology, rainfall, GHG emissions, trend maps on web mapping application (Table 3.6). - Bivariate Map Design: An effective bivariate map can communicate the relation between variables that is more essential than a set of individual maps of each attribute. Bivariate maps can enhance the ability of users to detect and comprehend important phenomena and major conclusions. This research selected important phenological parameters to show in bivariate maps. The important steps to design bivariate map are to: 1. Select appropriate variables/parameters to represent a relationship and the meaningful correlation, for example, the start date and the end date of the growing season. 2. Select appropriate techniques to show the pattern of phenological phenomena in each year and compare the change of these phenomena in different years, for example mean/median deviation with diverging scheme is an effective technique to show the increase and decrease pattern of vegetation condition (EVI) or demonstrate the changes in the start and end dates in the early or late growing season. This bivariate technique can answer the question where these patterns are located, and how much the magnitude of these attributes are. In addition, two 91 quantitative variables in this technique can show the correlation where both variables are high or low values, and where one variable is high while the other is low. 3. Select data classification that depends on data characteristics, such as using biweekly time-series data, the appropriate range should be one or two weeks to present the early and late growing season. Furthermore, no more than three classes are appropriate for bivariate maps because they can increase cognitive overload leading to users’ confusion with the maps. 4. Select color scheme to display the pattern of phenological information in each category. This research applied Brewer’s color scheme (1994) for bivariate mapping. For example, sequential/sequential scheme needs two hue and lightness differences to show high or low values, while diverging/diverging scheme bases on dark hues at each corner of the legend to represent categories that are extremes for both variables. A light and white color is placed at the center of the legend to represent a class that contains the critical value or midpoint of both variables. The remaining colors are lighter than the corners. - Animated Map Design: The prototype visualization system allows users to explore, monitor, and compare patterns and trends of vegetation phenology in SEA by using map animation. The major steps of map animation design are: 1. The first step is the same as bivariate maps to select appropriate factors: selecting appropriate parameters, data classification, the number of classes, and the color schemes, particularly simplifying the data before using in the map animation. The techniques used for animated maps are important to make animated maps effective and easy to read. 2. This prototype uses interactive temporal and spatial controls that are useful for focusing on details and facilitating an understanding of the process of change. When users can replay, pause, stop between scenes, users have more chances to review events they may have 92 missed within the sequence. Cox (1990) suggested that dynamic display were preferable for exploratory analysis of multivariate Earth science data. The difficulty in visual interpretation, the symbols, and the cognitive overload can be reduced by interactive animation. 3. Animated maps also require dynamic visual variables. Using dynamic visual variables in map animated design help to reduce the challenges of map animation such as change blindness by applying appropriate duration, rate of change, or order. 4. Other techniques can also enhance the effective animated map. For example, maps and text explanations appear close to each other and relate to one another in meaningful ways can reduce the amount of working memory available for learning. 5. The last factor is to carefully apply color on maps. It is particularly important in interactive and animated map contexts because the reader must attend to changing patterns on the maps and have little time to look back and forth to a map legend. 3.8.3 Visualization Implementation - Interactive Web Map: Interactive web map is used in this visualization system to create important functions and disseminate maps to users via Internet. Web mapping architecture consists of the following components (Figure 3.11): 93 Figure 3.11 Web mapping architecture 1. Client-side Component: Users can select parameters and display maps running in general web browsers through user interface (Mathiyalagan et al., 2005). 2. Server-side Component: Server side composes of three components: - Web Server: The function of the web server is to deliver web pages to clients by request. The web server communicates between clients and map server via gateway, for example, CGI-script (Gartner, 1999). 94 - Map Server: A map server provides the connection from the database to a web server. A map server serves mapping and spatial analysis functions as well as manages and distributes map data. One of the potential web map servers is University of Minnesota (UMN) MapServer. These server-side functions are sometimes referred to as web map services (WMS). The WMS is a standard protocol for serving geo-referenced map images over the Internet that are generated by a map server and a distribution service open to anybody needing information from the service, for example, Google map service. - Database Server: This server provides the Database Management System functionalities. It is responsible for searching the database and retrieving results from it. When users select parameters on a web browser, a web browser sends the request to a web server. Then, a web server interprets the request and passes it to a map server. The map server connects to database and processes the data. The results are sent back to the web server by converting into HTML format to display in the client-side web browsers (Mathiyalagan et al., 2005). Maps and images for the main web page in IPA were installed in the database and they are automatically and dynamically generated as users requested by MapServer. This database that contains images and maps were categorized in different folders based on the types of data and topics. In order to serve maps on the web, the UMN MapServer was used to create different map services for different data. Map services are determined by proposed functions (e.g., zoom, pan, parameter selection). One computer was selected as the web server to communicate between clients and the map server. The basic setup of a map server relies on two major components (Hachler, 2003). First, the Template File controls the display of a map server output in a web page (the design of the graphical user interface and how users can interact with 95 the application). Second, the Map File is the configuration file and controls all other aspects of a map server application such as the layers to display and the display parameters. - User Interface and Interactive Functions: The user interface of visualization systems interacts with users through a web browser. The user interface was developed using MapServer, HTML, Javascript, and python. The web user interface has menus to display, control, and link to information, e.g., phenological parameters as well as toolbar functions. There are two main web pages for this system. The first web page is the main page for this system, including all parameters and functions. The second web page is for map animation. The interactive functions enable users to change parameters and the application reacts on those changes. This system built two groups of interactive functions—basic and advanced functions—for users to explore maps and images. The basic functions provide navigation (zoom and pan tools), map coordinate, map backgrounds and map-based overlay, scale bar, and overview map using Google Maps API. Transparency adjustment was also included for users to adjust transparency of overlay thematic map. The advanced functions provide dynamic generating layer content, which update automatically corresponding to database. This function enables users to select data layers and years. Additionally, dynamic legends and time-series graphs were added to this system. For dynamic legends, when users select parameters such as EVI in 2001, the program processes the data and then creates the legend to show on a web browser. Therefore, the program automatically processes and shows the data as the developer defined. For graph creation, this system provides both static and dynamic graphs created by Google Chart API. The static graphs of EVI and rainfall as well as GHG emissions appear on one side of web page as referenced data. The 96 dynamic graph is created from the database and displays over the map. This tool allows users to click on any pixels on the map, and graph of that pixel will generate automatically. For example, when users click on one pixel of the start date of the growing season map, the system will generate graph of that pixel to show start date for ten years. This advanced function is very useful for users to explore and compare the patterns of that parameter over time. - Interactive Flash Map: Map animation of ten-year EVI and phenology were created in this visualization system. Maps and legends for map animation have to be created and prepared before building the animation. For map animation, the functions of map animation control, temporal legend, and temporal control were created by using Flash, Photoshop, and Illustrator. The series of completed maps in each phenological parameter were developed into animation and synchronized with the temporal legend and animation control by using Flash. Map animation in this prototype provides interactive animation and temporal control to represent changes over time and facilitate an understanding of the process of vegetation condition change in SEA. The tools in this prototype allow users to select the data, such as phenological parameters, to explore, monitor, and compare the patterns and changes of phenology. In addition, the start and stop time as well as duration of animation were set on the system. As a result, users can control animation by selecting play, pause, stop, as well as selecting any parameters on the tool bar to explore animation. 97 CHAPTER 4 SPATIO-TEMPORAL VARIATION OF VEGETATION PHENOLOGY IN SEA 4.1 Introduction The long-term change of phenology reflects the key aspects of global environmental impacts. This spatio-temporal extent of phenological change can be mapped with satellite timeseries data. The Indochina Peninsula of SEA was selected as the study area in this research. SEA is one of the world’s regions that is most vulnerable to the impacts of climate change and climate change has affected the ecosystem of this region. Therefore, this chapter discusses the spatial characteristics of the vegetation phenology in SEA by extracting EVI and seven phenological parameters and identified seasonal and inter-annual variability of vegetation phenology. The seven phenological parameters are the start (SGS), the end (EGS), the length (LGS), the middle, the amplitude of the growing season, the large integral, and the small integral. The patterns and trends of phenology were analyzed on a regional scale. Furthermore, specific hotspots were selected to illustrate the significant changes on a local scale for agricultural and naturally vegetated areas. In this research, the spatio-temporal variation of phenological parameters was investigated in both the first and the second growing seasons. EVI and phenological parameters were derived from MODIS products (16-day EVI with 250 m resolution from 2001to 2010). These satellite time-series images can provide an overall assessment of key ecological attribute changes over time across different ecosystems. Additionally, this chapter compares satellitederived phenology with field observations to assess the phenology results based remote sensing time-series data. 98 4.2 Seasonal Characteristics of Vegetation Profiles and Land Use Land Cover Types in SEA The seasonal profiles derived from MODIS EVI time series data exhibit distinct seasonal cycles with a well-defined greenness period that can be used to depict phenology and seasonal variations of vegetation activity (Huete et al., 2002). The rapid increase of the EVI profile in the beginning of the growing season and the rapid decrease of the EVI profile at the end of growing season make the timing, length, amplitude, and integral of the growing season easy to define. As a result, these EVI profiles have been widely used to delineate plant growing season. The seasonal EVI profiles were selected to show the growing season in order to understand phenological characteristics in SEA. The major land cover types in this region are forests and agricultural lands. Figure.4.1 shows the land cover in SEA in 2005. The tropical evergreen forests are found in the northern and southern regions of the Peninsula. The tropical deciduous forests are mostly located in the northern and eastern portions of the Peninsula. Land use in the whole region including Insular SEA has changed rapidly to agriculture in the past century (ADB, 2009). Agriculture is the most important sector in SEA, accounting for 11% of gross domestic product (GDP) in 2006 and providing 43.3% of employment in 2004 (ADB, 2009). The planted agricultural land areas are approximately 115 million ha in size. The important agricultural areas in this region are located in central and southern Myanmar, central and northeastern Thailand, northern and southern Vietnam, and central Cambodia. The main crops in this region are rice, maize, oil palm, natural rubber, and coconut. This region produced approximately 140 million tons of milled rice per year from 2002 to 2007 (ADB, 2009), has been a major supplier of grain, led by Vietnam and Thailand, and is the largest producer of palm oil and natural rubber in the world (FAO, 2007). 99 Irrigated Croplands Rainfed Croplands Croplands Deciduous Forest/ Grassland/Shrubland Evergreen Forest Urban Bare Areas Water Figure 4.1 Land cover in SEA from GLOBCOVER in 2005 100 There are multiple cropping frequencies (more than one growing season in a year) in this region, according to appropriate climatic conditions and plentiful water resources associated with the intensive farming management practices, and rice is the dominant crop in the second and third growing seasons. According to the research of Bridhikitti & Overcamp (2012), the rice paddy areas of SEA are approximately 30% of the world total. Rice paddy areas are mainly found in southern Myanmar, central and northeastern Thailand, and northern and southern Vietnam (Figure 4.2). There are four major rice ecosystems in this region: irrigated rice, rainfed lowland rice, deepwater rice, and upland rice. Figure.4.2 shows that the rainfed rice ecosystem is found in northeast Thailand, Central Myanmar, and Cambodia. The irrigated rice ecosystem is mostly observed in central Thailand, southern Myanmar, northern and southern Vietnam. The deepwater rice ecosystem is distributed in southern Vietnam. In addition, multiple rice cropping frequencies are also found in this region including single, double, triple (multiple) rice crops. Double and multiple rice crops are found in the major irrigated rice ecosystems (Figure 4.3). Seasonal EVI profiles of different land cover types are shown in Figure 4.4. These multitemporal VI profiles of a crop reflect the crop’s general phenological characteristics (Wardlow et al., 2007). The typical land cover types in SEA are forest, rice, and farm crops (e.g., cassava, corn, sugarcane). Evergreen forest shows very low amplitude due to a monotonous growing season (Xiao et al., 2009). In other words, there are no seasonal dynamics of leaf phenology in evergreen tropical forests. In contrast, deciduous forest has very high amplitude because of the high density of green cover and the amplitude of the deciduous forest is also higher than croplands. 101 Rice Cropping Ecosystems Deepwater Rice Upland Rice Rainfed Rice Irrigated Rice Figure 4.2 Rice cropping system in SEA 2007 Source: Bridhikitti & Overcamp (2012) 102 Rice Cropping Frequency Single Double Multi Figure 4.3 Rice cropping frequency in SEA 2007 Source: Bridhikitti & Overcamp (2012) 103 Furthermore, the number of growing seasons can be detected from the EVI profile to define one annual season or two annual seasons of vegetation in each pixel. The results show that most croplands have only one growing season throughout the year and the amplitude is slightly higher than that of rice. Also, different cropping frequencies were observed for paddy rice: single, double, and triple rice. The start and end of the growing seasons are also different for paddy rice from different locations. In addition to a clear signal of seasonal dynamics, the EVI profile can indicate the integral change and the timing shift when comparing with different years, and also represent the change of cropping frequency, for example, rice paddy changes from one to two growing seasons in a year (Figure 4.5). Crops phenology is likely dependent on agricultural practices, which may vary greatly among farmers as well as the variations of soil properties, water supply, and fertilizer management. Therefore, the EVI profile of each pixel has distinct seasonal dynamics indicating different phenological patterns (the timing, length, and amplitude) for different land cover types. Moreover, in the same pixel, the profile can be different from year to year due to climate variation or farming practices. Patterns and trends of EVI and phenology in the study area (Indochina Peninsula) on the regional and local scales are discussed in the following sections. 104 Figure 4.4 Seasonal EVI profile patterns 105 The number of season changes 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 2001 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 The 2007date and magnitude changes 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 LANDSAT 5 TM Vietnam 10/09/2005 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 LANDSAT 5 TM Thailand 11/21/2004 Figure 4.5 Seasonal EVI Profile Changes 106 forest rice rice 2001 2004 2007 4.3 Spatio-temporal Variation of EVI 4.3.1 Spatial Patterns of Interannual EVI The temporal EVI patterns in SEA from 2001 to 2010 were identified from the annual mean of EVI 16-day time-series composite images, which illustrated the inter-annual variability over the 10-year study period. The annual mean of the EVI is an indicator of the overall greenness of vegetation throughout the growing season. The results of spatio-temporal EVI variation depict the regional vegetation dynamics and long-term trend in light of climate and land cover changes. The spatial distribution of 10-year mean EVI from 2001 to 2010 (Figure 4.6a) exhibited larger values in highly productive agriculture areas and forests. These areas were mostly located in northern and southern portions of the Peninsula and in eastern Vietnam. Deciduous and evergreen forests were found in those areas. In contrast, regions where crop conditions were poor or occupied by low biomass and deforested areas tended to have lower EVI values. Relatively low EVI values were found in central Myanmar and northeastern Thailand. These lower EVI values are located in agricultural areas, compare to forested areas. The main agricultural areas in this region are in southern Myanmar, central Thailand, northern and southern Vietnam. The spatial distribution of EVI mean deviation from 2001 to 2010 is presented in Figure 4.7. Due to remarkable variation of EVI values in 2005 and 2010 according to El Niño event, mean EVI over ten years did not include 2005 and 2010. The results clearly indicate regional changes in the EVI values of each year from the 8-year mean EVI. Increases in EVI values were evident from 2001 to 2003 (wet year), followed by a marked decrease from 2004 to 2005 (dry year). The EVI increased gently in 2006 and 2007 and tended to decrease again to the 107 end of the analysis period in 2010. The years 2004-2005, and 2009-2010 were El Niño years with remarkably low rainfall, which apparently had a significant impact on the total vegetation growth in this region (TMD Thailand, 2010). The second increasing trend suggests the recovery of the ecosystems in the region after the extreme climate events. However, the second decreasing trend was found again at the end of study period, particularly in western Myanmar. These results indicate variable vegetation dynamics over SEA, mostly related to major climate events. 4.3.2 Trend Analysis The trend analysis of EVI was performed using the Mann-Kendall (MK) and Sen’s method to show the slope of the temporal trend for each pixel in this region (Figure 4.6b). The EVI profiles, which have the potential to represent seasonal and annual dynamics, show the total trends and changes of vegetation conditions for 10 years at a 90% confidence level (α=0.05). The positive and negative trends of EVI indicate the overall greenness of vegetation throughout the growing season, which are increased and decreased greenness, respectively. The results indicated that the overall regional trend of EVI in SEA was decreasing from 2001 to 2010. EVI patterns revealed a significant negative trend throughout most of SEA, particularly in western Myanmar, northern Thailand, and northern Vietnam. The areas which had a distinct positive trend of EVI were in southern Cambodia and Vietnam. In order to understand EVI changes in this region, some specific hotspots were selected to show the significant changes on a local scale. EVI trends of eight study sites, representing both agricultural areas and naturally vegetated areas, were investigated. The representative sites of forest areas are located on highlands of Mainland SEA including uplands of northern Myanmar, Thailand, Vietnam and Laos, and eastern Cambodia. These areas are subtropical forest eco-regions with sparse populations and low human disturbance (Tottrup et al., 108 2007). The major agricultural areas in this region, particularly rice paddy fields, are located on the river deltas and lower floodplains: the Irrawaddy River delta in southern Myanmar, the Chaopraya River delta in central Thailand, the Mekong River delta in southern Vietnam, and the Red River Delta in northern Vietnam (ADB, 2009; Xiao et al., 2006b; Bridhikitti & Overcamp, 2012). These areas are the major rice growing areas in the world (FAO, 2007; ADB, 2009). Therefore, the areas on the river delta are not only the major crop sources in this region but are also the largest sources of world food production. The eight typical locations with both naturally vegetated areas and agricultural areas exhibited the changes of EVI and land cover map (Figure 4.1) was also used to investigate these changes as follows: - Trends of naturally vegetated areas clearly showed decreased EVI in all selected sites (Figure 4.8). This pattern resulted from the decrease of forest areas in this region during the past ten years. Although some parts of this region still have a high forest cover such as northern Laos and Vietnam, some areas, the northern part of Thailand and the northern and eastern parts of the Peninsula, have relatively lower forest cover and the decreasing patterns of EVI were obviously seen in this site. This result suggests the forest loss in this region. - Most agricultural areas showed a decreasing trend of EVI (Figure 4.9). The large patches of negative trends were found in central Thailand and southern Vietnam. Specifically, rice paddies in central Thailand clearly exhibited the negative patterns. However, positive trends, which indicate the increase of vegetation greenness, were found in other croplands, e.g., cassava, and sugar cane in Thailand. In Vietnam, there were both positive and negative patterns. Most rice paddies in southern Vietnam have positive trends except for a single rice cropping system and rainfed rice, which showed a significantly decreasing trend. The 109 positive patterns were found in rice areas of northern Vietnam and southern Myanmar. Therefore, the directions of EVI trends in agricultural areas of this region were all spatially varied depending on farming management and climate variations. 4.4 Spatio-temporal Variation of Phenological Parameters 4.4.1 Spatial Patterns of Interannual Phenology The 10-year mean phenological values were calculated for each pixel in the MODIS EVI time-series data in SEA from 2001 to 2010 (Figure 4.10). However, pixels containing EVI profiles that show no or little annual variation (a lack of seasonal vegetation cycles) such as water, bare soil, tropical rain forests, were cut off from the analysis and shown in the white color. In this research, phenological parameters were calculated on an interannual basis, e.g., the data in January 2002 – July 2003 were calculated for the phenological results in 2002 (more details in Chapter 3 Section 3.4). The first growing season (Season 1) is defined as the first growing period of that year and the second growing season (Season 2) is referred to as the second growing period, which the start date of the second growing season is on the same year of the first growing season but the end date may be in the same year or in the next year. The spatial variations of seven phenological parameters are shown in Figure 4.10. 110 Figure 4.6 (a)10-year mean EVI (2001-2010) and (b) EVI trend (2001-2010) Figure 4.7 EVI mean deviation 2001-2010 111 EVI Trend North of the Peninsula East of the Peninsula % EVI / Year < -0.01 -0.01 - 0.0 0.0 0.0 - 0.01 > 0.01 Figure 4.8 EVI trend of hotspots in forest areas EVI Trend Central Thailand North of Vietnam South of Vietnam South of Myanmar % EVI / Year < -0.01 -0.01 - 0.0 0.0 0.0 - 0.01 > 0.01 Figure 4.9 EVI trend of hotspots in agricultural areas 112 Figure 4.10 10-year mean phenology, the first growing season (2001-2010) 113 Figure 4.10 (cont’d) - The First Growing Season (Season 1) 1. The start of the growing season (SGS): SGS usually ranged from March through May in SEA. SGS shifted toward the later dates from the east to west. The typical start 114 date from March to April was clearly shown in the northeast and the center of the Peninsula. However, the northwestern part of the Peninsula experienced slightly later starting dates in May. Most agricultural areas exhibited early start of the growing season, for example, rice paddy in northern Vietnam and central Thailand showed distinct early patterns. 2. The end of the growing season (EGS): The growing season ended at different times of the year depending on the type of the land cover and land use management practices, ranging from June to January of the next year. However, in general, EGS shifted toward the later dates from the east to west. In the eastern portion of the Peninsula, distinct areas were evident where the season ended earlier in June. The corresponding starting dates in these areas ranged from late January through March. This area is the large rice paddy area in northern Vietnam that has the same cropping pattern (has similar SGS and EGS) in the first growing season. Comparatively, the northwest portion of the Peninsula (Myanmar) experienced a later growing season, starting sometime in May and ending in January and February of the next year. The land uses of this area were the rainfed and irrigated croplands, and the timing of the growing season was different from the eastern portion of the Peninsula. 3. The length of the growing season (LGS): The average length of the growing season was approximately 180-240 days. LGS was generally longer in the northwest through the northern portion of the Peninsula. Conversely, a large part of the central and eastern portions experienced shorter growing season lengths. This trend was supported by the general pattern of starting and ending dates described above. The longer season lengths in the northwestern and northern regions are likely due to moisture-rich environments, forests, and perennial crops with longer growing seasons compared to the crops found in the east. 115 4. The middle of the growing season: The middle of the growing season usually ranges from May through September. Generally, the middle date of growing season in this region was in August. However, the middle date in northern Vietnam was earlier in May. The later mid date was found in central Myanmar. 5. The amplitude of the growing season: The mean amplitude, an indicator of the overall vegetation amount in the growing season, varied greatly across the Peninsula. The fairly low amplitudes were found in the agricultural and dry areas in the eastern and central parts of the Peninsula. The highest amplitudes were found in northern and northwestern portion of the Peninsula, where the vegetation development, evergreen and deciduous forests, are mixed. 6. The large integral: The large integral map, which represents high annual net primary production (annual biomass), shows high values in the forests of the northern and eastern Peninsula (eastern Myanmar, northern Thailand, Laos, and Vietnam, and eastern Cambodia), whereas low integrals were found in agricultural areas in central and southern Myanmar, central and northeastern Thailand, and northern and southern Vietnam. The results of this map are related to the amplitude of the growing season map. The areas that displayed high amplitude of the growing season also showed high values of total vegetation production. Although there are large forest areas in the southern portion of the Peninsula, the large integral in this area showed very low values. This is because there is no variation in seasonal vegetation cycles in this area due to the rainy tropical climate along with rain forests. Therefore, phenology cannot be extracted in this area. 7. The small integral: The small integral, which indicates seasonally active vegetation or seasonal change in net primary production (annual net growth), provided high 116 values in the same areas where the large integral was high. Some agricultural and dry areas, which received less rainfall, had lower values, suggesting low productivity. - The Second Growing Season (Season 2) The 10-year mean of each phenological parameter for the second growing season in SEA is shown in Figure 4.11. Distinct patterns of phenological parameters in the second growing season can be observed in the major irrigated croplands, dominated by paddy fields: southern Myanmar, central Thailand, northern and southern Vietnam. These four locations are the major agricultural areas in SEA and the largest source of world food production. 1. The start of the second growing season (SGS2): SGS2 was generally from July to December. Northern Vietnam exhibited a large patch of SGS2 in July and August, which was earlier than other locations. The SGS2 in southern Myanmar, central Thailand, and southern Vietnam generally appeared later in November and December. Some parts of southern Myanmar and central Thailand had SGS2 in July and August. 2. The end of the second growing season (EGS2): EGS2 ranged from December to April of the next year. EGS2 of northern Vietnam ended earlier in November and December corresponding with the earlier start date. 3. The length of the second growing season (LGS2): LGS2 was shorter than the first growing season and lasted approximately 100-180 days. This length is related to the length of rice growing season, which lasts an average of 120 days for the second growing season in this region (Office of Agricultural Economics, 2010) 4. The middle of the second growing season: The middle of the growing season was usually around October to January. However, northern Vietnam experienced the earlier midseason ranging from July to September. 117 5. The amplitude of the second growing season: The amplitude of the second growing season did not vary as much as the first growing season. The average of the amplitude was around 0.3-0.4 representing the rice ecosystem. 6. The large and small integral of the second growing season: The patterns of the large and small integrals were very similar in most of the study area. The reason is that the crops in the second growing season are mainly rice; as a result, most agricultural areas exhibited the similar patterns. 4.4.2 Trend Analysis The trends of phenology were processed by using the MK and Sen’s method. Trends were detected in seasonality of phenology for a period of 10 years from 2001 to 2010 at a 90% confidence level (α=0.05). Figure 4.12 and 4.15 show the significant overall trends of seven phenological parameters. The phenological trends did not clearly show when displaying the whole region. More detailed analyses were performed at specific sites and are discussed below. Positive and negative trends refer to later and earlier dates for the timing of phenology (e.g., SGS, EGS) and also refer to longer and shorter lengths of growing seasons as well as increased and decreased magnitudes for amplitude and integral. - The First Growing Season (Season 1) This section explains the changes and trends of phenological parameters: the start, end, length, and middle of the growing season, the amplitude, and the large and small integral by focusing on the first growing season. The overall phenological trends appeared to shift towards a later and slightly longer growing season in the Peninsula. Most of the study areas exhibited a later trend in SGS and EGS as well as a longer LGS (Figure 4.12). However, significant trends toward earlier SGS 118 were observed in western Myanmar, central, northeastern and eastern Thailand, and the eastern portion of the Peninsula. The overall changes of SGS, EGS, and LGS were within 2 weeks. The middle of the growing season also indicated the later trend. The amplitude and integral trends, which represent biomass production, illustrated decreasing trends in the eastern portion of the Peninsula (Figure 4.12). The trends of the large and small integrals showed a correspondence with the trends of EVI. The decreasing integral trends were found in western Myanmar, northern Thailand, northern Vietnam, and the eastern portion of the Peninsula. The areas which had a distinct positive trend of the amplitude and large integral were in southwest Myanmar and northeast Thailand. A more detailed analysis of phenological changes in this region was explored in eight representative sites (the same sites used in the EVI trend analysis section). The selected forest areas are the valuable forest areas with low level of human disturbance in this region (Tottrup et al., 2007) and the representative agricultural areas are the major sources of grain in the world (FAO, 2007, ADB, 2009). Therefore, the patterns and trends of phenology in these sites, which indicate the change of the growing season, would be useful for forest protection and agricultural management. Spatial patterns of phenological changes in naturally vegetated areas and agricultural areas were examined in eight locations and compared with land cover map (Figure 4.1) and rice ecosystem and rice cropping frequency maps in Figure 4.2 and 4.3 (Bridhikitti & Overcamp, 2012) and rice cropping pattern (Sakamoto et al.,2006): - In naturally vegetated areas, forest areas were found in regions that displayed both later and earlier trends. Forest areas in the western and eastern portions of the Peninsula exhibited an earlier trend, whereas forest areas in the northern portions showed a later trend 119 (Figure 4.13). Significant trends toward later SGS and EGS presented in northern Myanmar and northern Thailand, whereas significant trends toward earlier SGS and EGS were found in southern Laos and eastern Vietnam. Therefore, this eastern portion of the Peninsula is distinctly different from other forest areas in this region. This site showed earlier dates of growing seasons with decreased trends of amplitude and integral. Although other forest areas (i.e., the north of Myanmar and Thailand) had decreasing trends of amplitude and integral, the growing seasons changed to later dates. However, an increase of amplitude and integral with later growing season dates were found in western Myanmar. Moreover, there were also greater increases of the small integral, which indicates the increase of seasonally active vegetation or annual net growth at these sites. The change of timing of the growing season may be due to climate variation, while the change of amplitude and integral with the variation of the dates of the growing season in forest areas suggests land use change. This is because the production or magnitude not only changes but the timing also changes. In particular, if the small integral is also increasing, it is possible that forest has changed to cropland which has occurred in southern Laos and eastern Vietnam. - Most of the agricultural areas in this region showed later trends in SGS and EGS from 2001 to 2010 (Figure 4.14a). In central Thailand, irrigated rice exhibited a significant positive trend indicating the growing season started and ended later more than 2 weeks per year with a longer growing season. Rainfed rice demonstrated opposite patterns with earlier trends in start and end dates of the growing season and slightly shorter length. Other crops (e.g., cassava, sugarcanes) in Thailand started later and ended earlier with shorter lengths. The later pattern was also found in northern Vietnam. However, a significant negative trend presented in most of 120 southern Vietnam, which showed the growing season started and ended earlier and had a shorter length for irrigated rice. Interestingly, the rainfed rice in the southernmost part of Vietnam exhibited the later patterns for SGS, EGS, and longer patterns for LGS. The positive and negative trends were found in agricultural areas of southern Myanmar. Rainfed rice is likely to be related to rainfall but irrigated rice is more influenced by human management. The middle of the growing season was in agreement with SGS, EGS, and LGS trends. The amplitude and integral showed increasing trends in most agricultural areas, particularly in irrigated rice in central Thailand and rainfed rice in southern Vietnam. This result suggests the increase in crop production. - The Second Growing Season This section explains the changes and trends of phenological parameters in the second growing season: the start, end, length, middle of growing season, amplitude, large, and small integral. The trends of the second growing season were clearly presented in irrigated croplands located in four locations: southern Myanmar, central Thailand, northern and southern Vietnam. These four locations are important sources of world food production. Generally, the phenological trends of the second growing season in SEA tended to shift toward later dates with longer growing season (Figure 4.14b ande 4.15). SGS and EGS of the second growing season in central Thailand and southern Vietnam had become later with a longer growing season. The amplitude and the integral had increased. In contrast, the growing season in northern Vietnam, started and ended earlier with shorter length. The amplitude and integral had a tendency of decreasing. Similar to central Thailand and southern Vietnam, the results illustrate that rice paddies in southern Myanmar had a trend of later SGS and EGS, however, the amplitude and integral were decreasing. 121 Figure 4.11 10-year mean phenology, the second growing season (2001-2010) 122 Figure 4.11 (cont’d) 123 Figure 4.12 Phenology trend, the first growing season (2001-2010) 124 Figure 4.12 (cont’d) 125 Phenological Trend in the North of the Peninsula Start date Week / Year End date < -2 -2 - 0 0 0- 2 >2 Week / Year < -2 -2 - 0 0 0- 2 >2 Amplitude Length Week / Year EVI / Year < -0.03 -0.03 - 0 0 0 - 0.03 > 0.03 < -2 -2 - 0 0 0- 2 >2 Small Integral Large Integral Accumulated EVI / Year Accumulated EVI / Year < -0.3 -0.3 - 0 0 0 - 0.3 > 0.3 < -0.6 -0.6 - 0 0 0 - 0.6 > 0.6 Figure 4.13 Phenology trend of hotspots in forest areas 126 Start Date End Date Central Thailand North of Vietnam South of Myanmar South of Vietnam South of Myanmar North of Vietnam Central Thailand South of Vietnam Week / Year < -2 -2 - 0 0 0- 2 >2 Length Central Thailand South of Myanmar Amplitude Central Thailand South of Vietnam Week / Year North of Vietnam South of Myanmar North of Vietnam South of Vietnam EVI / Year < -0.03 -0.03 - 0 0 0 - 0.03 > 0.03 < -2 -2 - 0 0 0- 2 >2 Figure 4.14a Phenology trend of hotspots in agricultural areas (the first growing season) 127 Figure 4.14a (cont’d) Small Integral Central Thailand North of Vietnam Accumulated EVI / Year South of Myanmar South of Vietnam 128 < -0.3 -0.3 - 0 0 0 - 0.3 > 0.3 Start Date North of Vietnam Central Thailand End Date Central Thailand North of Vietnam South of Myanmar South of Vietnam South of Myanmar South of Vietnam Week / Year < -2 -2 - 0 0 0- 2 >2 Amplitude Length Central Thailand North of Vietnam Central Thailand North of Vietnam South of Myanmar South of Vietnam South of Myanmar South of Vietnam Week / Year EVI / Year < -0.03 -0.03 - 0 0 0 - 0.03 > 0.03 < -2 -2 - 0 0 0- 2 >2 Figure 4.14b Phenology trend of hotspots in agricultural areas (the second growing season) 129 Figure 4.14b (cont’d) Small Integral Central Thailand North of Vietnam Accumulated EVI / Year South of Myanmar South of Vietnam 130 < -0.2 -0.2 - 0 0 0 - 0.2 > 0.2 Figure 4.15 Phenology trends, the second growing season (2001-2010) 131 Figure 4.15 (cont’d) 132 4.4.3 Trend from Bivariate Mapping In addition to the trend analysis processed by the MK trend test, this research also applied the bivariate mapping technique to display two related phenological parameters on one map. The purpose of bivariate mapping is to compare and examine the relation between two attributes in order to represent individual distributions or the correlation between them. Three types of bivariate mapping were produced based on the related phenological parameters: start and end dates, length and amplitude, and length and large integral (Figure 4.16, 4.17, and 4.18). Each map was produced by calculating the difference of phenological values of each year with the 10-year mean values on a pixel by pixel basis. It is important to note that the legends of bivariate mapping represent important geographic phenomena. Specifically, the extreme values of both variables on the map can be seen in the four corners of the legend representing the patterns of change. The four important types of changes of the start and end date is defined as the following: 1. An early shift of the growing season (started and ended earlier) 2. An expansion of the growing season (started earlier and ended later, longer growing season) 3. A delayed shift of the growing season (started and ended later) 4. A shrink of the growing season (started later and ended earlier, shorter growing season) The four important types of changes of the length and amplitude, and length and large integral were defined as the following: 1. A shorter length with a decreased integral/amplitude 133 2. A shorter length with increased integral/amplitude 3. A longer length with increased integral/amplitude 4. A longer length with decreased integral/amplitude The results indicated distinct regional changes of two phenological parameters and can lead to the following major conclusions: - The Start and End of the Growing Season The spatio-temporal changes of the start and end of the growing season varied year by year (Figure 4.16). However, some explicit patterns can be discerned and general conclusions can be made. The growing seasons were mostly longer and shifted slightly to earlier dates in 2001. Northern Thailand and Laos experienced longer growing seasons during this year. From 2002 to 2004, the start and end dates shifted earlier with a shorter growing season length. The growing season in central Myanmar did not change much from 2001 to 2004; however the timing of the growing season exhibited the substantial shift to later dates with a longer length from 2006 to 2008. Northeast of Thailand showed a later shift of the growing season from 2006 to 2008 but the eastern portion of the Peninsular presented an earlier shift of the growing season in 2009 and 2010. The growing seasons in 2005 and 2010 were obviously different than those in other years. The growing season significantly shifted and changed in those two years; the growing season shifted to the later dates and has a considerably shorter length. The observation in 2005 and 2010 reflects the dry years, which were El Niño years in this region with remarkably low rainfall, particularly in 2009. These earlier and later growing seasons appeared to be in good agreement with the rainfall changes (see more details in section 4.2). These results suggest the relationship between the extreme climate events and growing seasons in this region. 134 - The Amplitude and Length of the Growing Season Visual interpretation of the amplitude and length of the growing season in the bivariate maps reveal that the general patterns over the period of ten years in this region showed decreasing amplitude (Figure 4.17). The length of growing season was clearly shorter in the first five years (2001-2005) and longer in the last five years (2006-2010). The amplitude exhibited significant decreasing patterns in 2006, 2008, and 2009. However, an increase in amplitude was found in 2004, 2007, and 2010 with a shorter growing season in Thailand and longer growing season in Myanmar. The increasing amplitude with a shorter growing season length was more pronounced in 2010, especially in central Thailand and Myanmar, but most of Vietnam experienced the decreasing amplitude with shorter growing season lengths during this year. - The Large Integral and Length of the Growing Season Most of SEA showed two distinct patterns of the large integral and the length of growing season during the 10-year period (Figure 4.18). A decrease in the large integral with a shorter growing season length was found in 2002 to 2005 and 2010. Conversely, an increasing large integral with longer growing season lengths were exhibited in 2001 and 2006 to 2009. In addition, forest areas had an increasing integral and longer growing season lengths in 2001. This pattern changed to decreasing integrals and shorter lengths from 2002 to 2007 and switched back again to the higher integral and longer length in 2008, particularly in northern Thailand, Laos, and Vietnam. The overall patterns in 2010 clearly show decreasing integrals with shorter growing season lengths. The changes of these phenological parameters presented on bivariate maps generally suggest a significant impact of severe climate on the total vegetation growth and the growing season, which was clearly observed in the growing season that shifted to the later dates 135 and has a considerably shorter length in extreme climate years (2005 and 2010). Moreover, the amplitude and the large integral also showed decreasing patterns in 2005 and 2010. The amplitude and the large integral are related to the annual biomass and net primary production. The changes of these phenological parameters reflect the higher or lower biomass in each year. When considering specific sites for the changes of amplitude and large integral together with the length of growing season, the results can not only indicate the changes of biomass but they are also possible to suggest the land use and land cover change in those areas. For example, in northeast Thailand, the croplands showed an increased amplitude and large integral combined with a longer growing season length in the last five years of the study period (2006-2010). These results indicate land use changes in this location. The agricultural report in Thailand (Office of Agricultural Economics, 2010) supports this conclusion: most of this area has changed crop type, e.g., changed from rice to cassava. Therefore, these phenological changes provide the understanding of ecosystem dynamics from regional to local scales. 136 Figure 4.16 Bivariate map: mean different between the start and the end dates (2001-2010) 137 Figure 4.17 Bivariate map: mean different between the length and the amplitude (2001-2010) 138 Figure 4.18 Bivariate map: mean different between the length and the large intergral (2001-2010) 139 4.5 Comparison between MODIS Phenology and Field Observations This research compared the phenological parameters derived from MODIS EVI with field data in order to observe errors and quantify the agreement between satellite derived phenology and field observations. The field data (two hundred points) were obtained by interviewing famers and head of villagers in central and northeast Thailand. The results of the comparison are shown in Figure 4.19. The SGS and EGS of the first 2 growing season plots suggested high correlations of SGS and EGS in Thailand with R values of 0.6 and 0.9, respectively. The SGS and EGS of the second growing season derived from MODIS 2 also showed strong agreement with field data (R = 0.9 for both SGS and EGS). However, we can notice that SGS and EGS obtained from MODIS tended to occur later than SGS and EGS from field data. In addition, the MODIS date varied but the date from field observations appeared unchanged. There is an expected difference between MODIS phenology and field data due to the differences between observation techniques. MODIS phenology is extracted from integrated canopy greenness, while the field data were collected from interviewing farmers and in some cases farmers provided only the estimated dates of the growing season, for example the st beginning of the month was assumed as the 1 of the month. However, this research considers these data to be sufficient for the comparison of the regional phenological patterns because: 1) the field data were obtained from the major crops in this region and 2) the MODIS pixels were selected that contained only one crop type or from large areas of only one crop type. The results indicate that MODIS phenology provides reasonable results and can be used to estimate regional phenological patterns in SEA. 140 Start Date of the First Growing Season 176 points Start Date of the Second Growing Season 200 400 R² = 0.669 R² = 0.899 350 150 MODIS (Doy) MODIS (Doy) 67 points 100 50 300 250 200 150 0 100 0 50 100 Field (Doy) 150 200 100 End Date of the First Growing Season 350 400 500 R² = 0.892 400 200 250 300 Field (Doy) End Date of the Second Growing Season R² = 0.857 MODIS (Doy) MODIS (Doy) 500 150 300 200 100 400 300 200 100 0 0 100 200 300 Field (Doy) 400 100 150 200 250 300 350 400 450 500 Field (Doy) 500 Figure 4.19 Comparison between MODIS phenology and field data in Thailand 141 4.6 Discussion and Conclusions 4.6.1 Spatio-temporal Variation of Vegetation Phenology in SEA The results presented in this chapter demonstrate that the annual cycle of vegetation phenology inferred from remote sensing can identify spatio-temporal patterns of the growing season at annual time scales. In terms of spatial difference, phenological patterns spatially varied depending on geographic locations and ecosystems as well as farming management and climate variations. The homogeneous patterns of phenology were found in the same ecoregion. The phenological dates (e.g., SGS, EGS) of natural vegetation appeared later in natural vegetation but started earlier in agricultural areas and forest areas had longer growing seasons than agricultural areas. Due to multiple cropping frequencies in this region and rice is the dominant crop in multiple growing seasons, this research presented the patterns of phenological changes in both the first and second growing seasons. Rice fields in a double cropping system with irrigation practice started earlier in the beginning of the years because most of the rice fields in this area have more than one growing season in a year, resulting in an earlier start in the first growing season. This pattern of earlier start was obviously seen in central of Thailand, northern and southern of Vietnam. In contrast, other croplands (e.g., cassavas, sugarcane) started their growing season in approximately May or June, which was later than the rice fields. This is because these crops need the rainfall that begins around May in this region for vegetation growth. In addition, the same crop types are spatially different in each location and may have differences in phenological characteristics due to different soil, water, climate, management, etc. Therefore, the timing of the start and end date of the growing season, and the 142 length of the growing season, as well as the magnitude of greenness (e.g., maximum VI) vary from place to place according to environment and management of each location. Rice fields are one of examples, which have different phenological characteristics in each location. However, different characteristics can be observed from irrigated and rainfed rice fields. Phenology is not only different in different land uses or locations; it can also indicate temporal change or shift due to climate variability or human management. In case of temporal analysis, the results indicate that the EVI in SEA from 2001 to 2010 was decreasing. Phenological trends exhibited a later and slightly longer growing season. The biomass production (the amplitude and integral) demonstrate decreasing trends in the eastern portion of the Peninsula. Regional vegetation dynamics over SEA illustrated a relative correspondence with major climate events such as El Niño in 2005. Various parameters showed distinctly different patterns in years of these extreme climate events. To obtain more information about EVI and phenological changes in this region, the representatives of forest areas and agricultural areas in this region were selected to investigate the trend and change of growing season. EVI trends representing overall greenness of vegetation showed the forest loss in the northern part of Thailand, northern and eastern portions of the Peninsula. Furthermore, the overall forest areas exhibited the decreases in amplitude and integral indicating the loss of high-quality forests, The major causes of forest loss are the shifting cultivation (small agricultural fields) and commercial logging, as well as large-scale infrastructure developments (Tottrup et al, 2007). In addition, the extreme climate events have negatively affected the forest, e.g., forest fires in the dry year (ADB, 2009). Therefore, the results of EVI and phenological changes in forest areas demonstrate the influence of human disturbance and climate variation. In agricultural areas, rice paddies in central Thailand showed decreasing 143 trends of EVI but increasing trends were found in Vietnam and Myanmar. The amplitude and integral of irrigated rice in central Thailand and rainfed rice in southern Vietnam indicated the increase in annual net primary production. The amplitude and integral of agricultural areas had an increasing trend suggesting the improvement of agricultural practices. Therefore, forest areas in this region tended to change according to climate variation and human disturbance, while trends of agricultural areas varied depending on farming management. The drivers of phenological change in this region are investigated in the next chapter. Information of vegetation dynamics derived from remotely sensed data is essential for local and regional natural resource managements. Understanding these seasonal patterns of vegetation activity is valuable for the characterization of vegetation, studying the impact of climate change and inter-annual variability, monitoring land degradation, and detection changes in land use and land cover. In addition, phenological shifts and changes could result in large changes in annual gross primary production and affect the carbon cycle, water cycle and energy fluxes through photosynthesis and evapotranspiration (Xiao et al., 2009). These effects consequently may have influenced on food security, water resources availability, and climate. More importantly, changes in average climate conditions and climatic variability will have a significant effect on phenology again as the feedback loop. As a result, phenological detection based on remote sensing can provide phenological information in space and time that is benefit the study of ecological terrestrial system and environmental management. 4.6.2 Challenges and Limitations Although EVI time-series data at an intermediate spatial resolution and high temporal resolution provide useful land surface phenology for monitoring and assessing spatial and temporal variations in vegetation amount and condition, there are a number of challenges 144 that still need to be addressed. It is important to note that the data sets and the methods described in this research are the first attempt to extract phenology in this region on a regional scale; consequently, the following issues should be considered when applying EVI time-series data for phenological detection in this region. - Data Quality MODIS EVI can provide enhanced temporal detection of land surface phenology but this product is still constrained by cloud contamination and atmospheric effects. The 16-day compositing period of the MODIS EVI product series collection 5 (MOD13Q1) are improved by using a new quality based filtering scheme for cloud and aerosol correction, a modified compositing method to deal with residual and mislabeled clouds, as well as applying EVI2 as the backup algorithm for cloudy pixels (Didan & Huete, 2006). The compositing method is processed by selecting the date providing the smallest view angle with the highest VI values for minimizing the view angle effects with the best atmospheric conditions (Didan & Huete, 2006). This product is expected to have a tremendously positive affect on post-processing. However, vegetation in tropical climates is most active during the rainy (i.e., cloudy) season resulting in the selection of clear pixels, which is greatly impeded by residual clouds in the compositing method. It is important to note that optical remote sensing in a moist tropical region with monsoonal climate faces challenging issues such as frequent cloud cover in the wet season (growing season) and fire-induced aerosols in the dry season. To deal with these problems, MODIS EVI products used in this research have applied the Savitzky-Golay smoothing function in the TIMESAT program in order to replace outliers, spikes, or missing values. However, this approach is sometimes insufficient to eliminate cloud contamination, particularly to detect phenology in tropical climates. This leads to the 145 unsuccessful extraction of phenological information for some pixels. Because of the concerns about data quality issues, annual data quality assessment maps were produced in order to compare the level of quality associated with annual EVI and annual phenological patterns (Appendix C). The data quality maps were processed by using reliability layers provided in the MOD13Q1 product. It is found that the lower data quality is enhanced during the growing season (May-Aug) due to the monsoonal climate. Therefore, it is impossible to avoid low quality data in the regional scale analysis. For the local scale analysis, this study attempted to select locations having high data quality. However, some locations with low data quality were also selected because they are the representative locations of this region. The interpretation of phenology in the areas of marginal quality data should be an important concern. This research tried to examine the phenological results in locations with marginal quality data. One of interesting location is in the Red River delta of northern Vietnam. This location showed low quality data in most of images in a year (There are 23 images per year for the MODIS EVI composited product.). The results of this study found two growing seasons of rice in this location and other studies also confirmed that most of rice in the Red River delta is irrigated with two crops per year (Young et al., 2002; Cooke & Toda, 2008; Vietnam Trade Promotion Agency, 2008; Bridhikitti & Overcamp, 2012). Therefore, it seems that even with low quality data, satellite derived phenology can provide reliable results of general phenological patterns; however, timing of growing season may be varied due to data quality. MODIS AQUA product is another alternative for the study of vegetation condition. However, data quality and biophysical conditions of the study area are important for the selection of data products and sensors in phenological study. 146 - Validation Validation is a key issue in remote sensing-based studies of phenology on a regional scale (Zhang et al., 2003). The best way to evaluate results is to use field data. While a variety of field programs for monitoring phenology have been initiated, these programs provide data that are typically species-specific and are collected on scales that are not compatible with coarse resolution remote sensing (Reed et al., 2009). In addition, ground phenological databases have not yet been satisfactorily validated due to the difficulty in obtaining sufficiently extensive ground observations throughout the growing season. This research has employed the comparison between phenological parameters derived from MODIS EVI with field data in order to observe errors and quantify the agreement between satellite derived phenology and field observation. There is no availability of phenological stations or ground-based phenology observation in this region. Therefore, the field data in this research were recorded by interviewing famers and head of villagers in central and northeast Thailand. As a result, field observations cover small areas in Thailand due to time and financial constraints. Although the results of the comparison showed good agreement between field observations and satellite estimations, the technique caused deviations between the phenological estimates and the ground data in some cases. For example, estimates of growing season dates from MODIS were later than those from field data and the MODIS dates varied but the dates from the field observations appeared unchanged. It is noteworthy that substantially more field data are required to fully assess the phenology based remote sensing time-series data. Additionally, ground surveys of phenology in a variety of environments are required and field data collection should be specifically designed to validate remote sensing phenology. This means that the field approaches should cover large area estimates of phenology 147 over heterogeneous study areas, rather than plant specific measures. In addition, as Reed et al. (2009) suggested, the dates of phenological record should be related to continuous vegetation indices (e.g., start and end dates of the growing season). Such an integrated effort would be logistically difficult and expensive, but would be extremely beneficial to remote sensing-derived phenological applications as well as global change studies. The second solution is to develop and expand ground-based phenology observation networks, such as the USA National Phenology Network (Reed et al., 2009). Another possible way to validate phenology based remote-sensing time-series data is to compare the results with high resolution images, for example, comparing MODIS time-series data with LANDSAT data. However, this technique is difficult due to the need to acquire all time-series LANDSAT imagery for one year and it requires the appropriate method for data filtering and phenological extraction. Furthermore, for regional scale analysis, it is considered to be time consuming and requires high computing performance. - Resolution MODIS satellite time-series data provide high temporal resolution and is free of charge; however, the 250 m resolution leads to the spatial heterogeneity of land cover classes in one pixel. In the large MODIS pixel over diverse vegetation types, each pixel reflects the integrated response across landscapes with diverse species and phenological behavior. These mixed pixels introduce error into the analysis of interannual variability. The interpretation of remotely sensed phenological parameters at this spatial resolution should be addressed. The high resolution satellite images (e.g., LANDSAT) offer an appropriate local scale analysis but there are some limitations with these data for effective phenological extraction. Although there are a number of limitations and uncertainties to estimate satellite phenology derived from MODIS time-series data, the method used in this research provides 148 important phenological information; in particular, patterns and changes can be observed to understand regional vegetation variability and ecosystem dynamics in this region. 149 CHAPTER 5 UNDERSTANDING THE DRIVING FORCES OF PHENOLOGICAL CHANGES 5.1 Introduction The results of chapter 4 indicated the change of vegetation phenology in SEA. Moreover, the EVI and phenological patterns varied spatially. Regional vegetation dynamics in this region exhibited patterns that were associated with major climate events such as El Niño. However, human modification is a possible driver of such changes. The objective of this chapter is to analyze the causes/drivers of phenological changes from both climatic and anthropogenic factors. The first part of this chapter focuses on the relationship between phenology and climate variability to understand environmental changes in this region because this relationship is a key to identify the influence of climate change on biophysical characteristics. Since ecosystem of SEA is highly sensitive to rainfall variation and most of SEA suffers from recurrent droughts and floods, which strongly affect their agricultural areas, the relationship between rainfall variation and phenological changes was examined in this chapter. Spatio-temporal patterns and trends of rainfall seasonality were analyzed by using TRMM rainfall (TRMM-3B42-daily product from 2001to 2010). Then, the response of vegetation phenology to rainfall seasonality and rate of phenological changes, with respect to rainfall seasonality in this region, were explored. This study will be helpful in understanding the climate in SEA and how it affects the vegetation. This knowledge is important for future study of global climate change. Additionally, land use data and rainfall variation on a local scale analysis were also used to indicate possible drivers of change between anthropogenically driven land cover change and interannual climatic fluctuations. The analyses on both regional and local scales are essential for understanding the effects of climate variability and human management on ecosystem dynamics in this region. This chapter also 150 illustrates the agreement between TRMM rainfall and station rainfall. The last section of this chapter discusses the influence of climate change and human management on phenological patterns in SEA. 5.2 Response of Seasonal Vegetation Dynamics to Climate Variation in SEA 5.2.1 Spatio-temporal Variation of Rainfall Seasonality in SEA Spatial distributions of rainfall were generated to represent 10-year mean annual rainfall and rainfall seasonality (the start of rainy season (SRS), end of rainy season (ERS), and length of rainy season (LRS) by using TRMM satellite time series data with the biological criteria (Jutakorn, 2011) from 2001 to 2010 (Figure 5.1). A high intensity of rainfall was found in the southern Peninsula and southwestern Cambodia. In contrast, low rainfall was observed in the northwestern Peninsula. The average amount of rainfall in most of SEA was approximately 1,000-2,000 mm/year. The rainy season in SEA generally started in April or May. The southern portion of the Peninsula tended to have more rainfall and an earlier SRS due to year round rainfall resulting from its rainy tropical climate. The average ERS occurred from September to December, with an average LRS of 180-240 days with the exception of the southern peninsula, which exhibited a later ERS and longer rainy season. The trend of rainfall seasonality was calculated by MK and Sen’s method. The interannual change results show only pixels with a significant trend (Figure 5.2). The northwestern portion of the Peninsula and most areas of Laos presented a clear negative trend in rainfall, indicating decreasing rainfall in this region. Strong negative trends were found in the western Peninsula and in southern Laos. The rainy season in this region generally tended to start earlier (0-7 day/year), end later (0-7 day/year) and had a longer length (0-7 day/year), particularly in Thailand. Distinctive rainy season patterns were evident in Myanmar, as the 151 season tended to start later and had a shorter duration. Conversely, the southern portion of the Peninsula experienced an earlier start of rainfall and longer rainy season. 5.2.2 Comparison between Mean Annual EVI and Mean Annual Rainfall The temporal mean EVI and mean annual rainfall show decreasing patterns in SEA from 2001 to 2010 (Figure 5.3). Increases in mean EVI were evident from 2002 to 2003, followed by a marked decrease from 2003 to 2005. The EVI increased again in 2006, then gradually decreased and continued in this trend through the end of the analysis period in 2010. Mean annual rainfall showed a decreasing trend from 2001 to 2010, with a significant decreasing trend in 2005 and 2010, which corresponded to the trend of mean EVI. The years of 2005 and 2010 were El Niño years with remarkably low rainfall, which apparently had significant impacts on total vegetation growth in this region. The increasing trend of EVI and rainfall was basically reset after the El Niño event. The second increasing trend suggests the recovery of the ecosystems in the region after the extreme climate event in 2003 and 2006. However, this region exhibited continuously decreasing EVI and rainfall patterns from 2007 to 2010. EVI and rainfall in this region exhibited similar patterns indicating the influence of rainfall variability on total vegetation growth. Therefore, variable vegetation dynamics over SEA are more likely to be related to major climate events. 5.2.3 Relationship between Seasonal Rainfall Fluctuations and Phenological Parameters Phenological patterns derived from MODIS products and rainy season observations from TRMM data over the ten year study period were used for investigating the response of vegetation to rainfall and identifying the sensitivity of various vegetation formations 152 to climate variability. The relationship between important parameters was examined: SGS-SRS, EGS-ERS, LGS-LRS. EVI and rainfall, and SGS and SRS showed a strong positive correlation in most of the study area (Figure 5.4). EVI and rainfall exhibited a positive correlation in northern, northwestern, and eastern portions of the Peninsula. SGS and SRS were highly correlated in Myanmar and some parts of Thailand, Laos, and Vietnam. EGS and ERS showed strong positive trends in Myanmar and Thailand. A significant positive trend in LGS and LRS were also apparent in some parts of Myanmar, Thailand and Vietnam. These results indicate that positive correlations were found in water deficient ecosystems, where ecosystems are dependent on rainfall for production. In contrast, a negative correlation was found in scattered pixels, where irrigated agricultural areas exist. 5.2.4 Rate of Phenological Changes with Respect to Rainfall Seasonality Trends of relation presented the rate of shift were performed by considering phenology with respect to the change in rainfall during 2001-2010 (how many days phenology changes when rainfall increases or decreases). The regression analysis was processed for pixels with significant correlation and with a corresponding increase or decrease of rainfall parameters for a period of 10 years. Phenological parameters with respect to corresponding change in rainfall parameters show obvious results in some cases and only few pixels showed significant trends of relation. Therefore, this research selected three important cases to show the rate of change for SGS-SRS, EGS-ERS, LGS-LRS. The spatial patterns can be detected as shown in (Figure 5.5). 153 Figure 5.1 10- year mean rainfall (2001-2010) 154 Figure 5.2 Rainfall trends (2001-2010) 155 Figure 5.3 Mean annual EVI and mean annual rainfall 156 Figure 5.4 Relationship between phenology and rainfall seasonality 157 Figure 5.5 Rate of phenological change with respect to rainfall seasonality 158 - SGS showed a positive trend with respect to SRS. This suggests that SGS shifted backward at rate of 1-3 days per one day decrease of SRS in the eastern portion of the Peninsula. - EGS showed a positive trend with respect to ERS. This indicates that EGS shifted forward at a rate of 1-3 days per one day increase of ERS in the northeastern Thailand. - LGS showed positive trend with respect to LRS. This demonstrates that LGS shifted forward at a rate of 1-2 days per one day increase of LRS in some parts of Thailand. 5.3 Driver of Phenological Change: Local Scale Analysis Changes and trends of vegetation phenology can provide strong scientific evidence for ecosystem dynamics. Both climate variations and human management are possible drivers of phenological changes. To understand phenological changes as explained in the previous section, the drivers of changes should be addressed. Therefore, this section analyzes the possible drivers of phenological changes between anthropogenically driven land cover change, and interannual climatic fluctuations on a local scale analysis (in Thailand). Although overall phenological change at a regional scale is important for examining drivers of change, local scale analysis can provide more precise and meaningful information about those drivers. In this section, the relationship between phenology and climate variability and land use/land cover changes were investigated in specific sites by selecting hotspots based on significant relationships among phenology and rainfall, striking trends, and differences of land cover characteristics. Four case studies were selected for this analysis: Site I and II in the northern Thailand and Site III and IV in central Thailand. First, the relationship between phenology and climate variability was explored to monitor and assess environmental changes in this region. This relationship is a key to understand the 159 influence of climate change on biophysical characteristics. As rainfall is the key factor to control vegetation growth duration in this region (Kramer et al., 2000; Cleland et al., 2007), this research focused on the response of phenology to rainfall variation in order to understand the sensitivity of various vegetation formations to rainfall changes. Phenological patterns were derived from MODIS products and rainy season observations were obtained from TRMM data over the ten year study period (as described in the previous section). Second, land use data were also considered to indicate the possible causes of change in phenology, which are human management or climate. Land use data for two periods in each site were provided by National Land Use Dataset of Thailand from Land Development Department, Ministry of Agriculture. This research also compared the phenological change and land use change with the rice ecosystem and rice cropping frequency maps (Figure 4.2 and 4.3). These maps were supported by Mahasarakam University in Thailand (Bridhikitti & Overcamp, 2012). 5.3.1 Site I and II in Northern Thailand The first two sites are located in northern Thailand. The first site is in Sukhothai province, which is mainly covered by forest, with small areas of rainfed rice and corn (Figure 5.6). The relationship between phenology and rainfall seasonality at this site showed a positive correlation in all pairs of parameters (SGS-SRS, EGS-ERS, and LGS-LRS). SGS displayed a later trend while EGS displayed an earlier trend; resulting in a shorter LGS for this site. Land use changes were also considered to identify the influences of human management on phenology change. Land use change and change detection were processed by using land use data in 2002 and 2009. The results show that land use was not changed substantially at this site (Table 5.1). The obvious change was that sugarcane and corn had increased, while perennial, and mix field 160 crop had decreased. Therefore, phenological changes at this site were mainly due to the changes in rainfall patterns. The second site is also located in northern Thailand, in the Uttaradit province. The land cover in this site is mainly rainfed rice with some large areas of evergreen and deciduous forests (Figure 5.6). Phenology and rainfall also showed positive correlations in SGS-SRS and LGS-LRS. When considering land use change in 2001 and 2009, the areas of cropland showed little change (Table 5.2). The change detection shows that a large portion of sugarcane was converted into cassava. The results of these two sites exhibited a strong relationship between phenology and rainfall and there was relatively little land use change. Consequently, change in rainfall is the major influencing factor for phenology at these two sites. Table 5.1 Land use change for site I Sukhothai province in 2002 and 2009 Land Use Urban Paddy Field Corn Sugar cane Cassava Mix field crop Perennial Orchard Evergreen Forest Deciduous Forest Water Miscellancous Land 2002 (%) 0.74 2.86 0.98 0.31 0.00 10.54 8.59 9.07 29.07 36.73 1.09 0.00 2009 (%) 1.60 3.29 8.44 0.74 0.03 0.31 3.30 9.45 28.35 42.85 0.99 0.65 161 % Change 114.27 14.72 760.95 136.61 2.52 -97.07 -61.54 4.20 -2.47 16.67 -9.53 65.23 Table 5.2 Land use change for site II Uttaradit province in 2001 and 2009 Land Use Urban Paddy Field Corn Sugar cane Cassava Mix field crop Perenial Orchard Horticulture Aquacultural Land Evergreen Forest Deciduous Forest Water Miscellancous Land 2001 (%) 11.11 38.40 3.97 4.94 1.58 8.39 4.27 2.15 0.00 0.00 11.84 7.79 2.91 2.65 2009 (%) 13.38 37.37 1.90 3.68 3.17 4.18 5.67 2.54 0.40 0.03 12.81 7.00 2.87 5.00 % Change 20.42 -2.70 -52.05 -25.37 100.93 -50.16 32.75 18.20 40.48 3.47 8.19 -10.17 -1.67 88.48 5.3.2 Site III and IV in Central Thailand Sites III and IV are located in central Thailand. The central area of Thailand contains large areas of paddy fields and a mixture of croplands. It is not surprising that most pixels in this area showed a negative correlation between phenology and rainfall seasonality because the croplands are highly affected by human management. Site III, located in Supanburi province contains mostly irrigated rice and some sugarcane (Figure 5.7). Phenology and rainfall showed a negative correlation here. This site is largely covered by irrigated croplands; consequently, it is affected by human management. In addition, change detection analysis of land use in 2000 and 2010 clearly indicated that the important land use change in this site was the areas of sugarcane, which were increasing and changed from paddy rice (Table 5.3). It can be concluded that the phenology in this site was very dependent on agricultural practices, which may vary greatly among farmers. Therefore, the 162 precipitation changes did not show a strong relationship with the phenological changes. On the basis of these observations, this site reflects human management rather than climate influences. The last site is also located in Lopburi province. Half of this site is covered by paddy fields and another half by sugarcane (Figure 5.7). Most of paddy fields in this site are rainfed rice. Although this site is in central Thailand and is mainly covered by croplands, phenology and rainfall were positively correlated. When exploring land use changes in 2001 and 2010 (Table 5.4), a noticeable change was the area of corn, which was converted to sugarcane. This site indicates that although land use is changed, the rainfall pattern also accounted for phenology changes, suggesting some croplands also need rainfall to control vegetation growth. In summary, natural vegetation is likely dependent on rainfall for vegetation growth duration. In addition, rainfed croplands also need rainfall for the growing season. Therefore, the results demonstrate the strong relationship between phenology and rainfall in these areas. However, most of the croplands, particularly irrigated crops, are more influenced by human management. Thus, phenology and rainfall exhibited a negative correlation. These characteristics need local scale analyses site-by-site to indicate the possible drivers of changes and to identify the hotspots or sensitive areas of ecosystem dynamics. 163 Table 5.3 Land use change for site III Suphanburi province in 2000 and 2010 Land Use Urban Paddy Field Sugar cane Cassava Mix field crop Perennial Orchard Horticulture Aquacultural Land Deciduous Forest Water Miscellancous Land 2000 (%) 6.09 84.73 3.92 0.00 0.22 0.00 1.85 0.00 0.00 1.16 1.24 0.79 2010 (%) 12.13 65.54 9.98 0.05 0.21 0.31 3.94 0.06 0.69 1.41 4.09 1.58 % Change 99.31 -22.65 154.20 4.94 -2.18 31.24 112.54 6.33 68.76 21.62 230.77 100.13 Table 5.4 Land use change for site IV Lopburi province in 2002 and 2010 Land Use Urban Paddy Field Corn Sugar cane Cassava Mix field crop Perennial Orchard Horticulture Aquacultural Land Deciduous Forest Water Miscellancous Land 2002 (%) 2010 (%) % Change 2.84 6.58 132.08 50.19 46.74 -6.88 39.15 0.34 -99.14 1.36 30.74 2163.07 0.00 3.77 377.03 0.00 2.41 241.11 0.32 2.73 761.28 4.71 2.36 -49.89 0.00 0.03 2.78 0.06 0.29 397.47 1.07 1.73 61.86 0.31 1.64 430.73 0.00 0.63 63.17 164 Figure 5.6 Land use changes in Site I and II in northern Thailand 165 Figure 5.6 (cont’d) 166 Figure 5.7 Land use changes in site III and IV in central Thailand 167 Figure 5.7 (cont’d) 168 5.4 Comparison between TRMM Rainfall and Station Rainfall This research compared TRMM rainfall with station rainfall. Monthly station rainfall and the average ten-year mean of rainfall seasonality (SRS, ERS, and LRS) from Thailand research (Jutakorn, 2011) were used for this comparison. As indicated from Figure 5.8, the results reveal that annual rainfall extracted from TRMM was strongly correlated with rainfall station data with 2 R values of 0.6-0.8 from 2001 to 2010. Furthermore, the results showed that the slope was lower than the 1:1 line, indicating that the TRMM estimation was lower than what was observed by rainfall stations, particularly for heavy rainfall. Previous research also mentioned to this underestimated rainfall patterns from TRMM and this problem is more pronounced for heavy rain rates (Huffman et al., 2007; Chokngamwong & Chiu, 2008). This issue should be addressed when applying satellite-precipitation sensors. For rainfall seasonality, Figure 5.9 clearly illustrated good correlation of SRS, ERS, and LRS between the rainfall parameters derived from TRMM and obtained from TMD during a 2 period of 10 years with R values of 0.5, 0.8, and 0.6, respectively. However, the discrepancy was also examined. The rainfall seasonality of TRMM showed earlier dates than what was estimated from TMD. The overall results indicate that rainfall from TRMM can be used to predict the rainfall seasonality. 169 Mean Annual Rainfall 5000 2001 5000 TRMM (mm) TRMM (mm) 6000 4000 3000 2000 R² = 0.807 1000 3000 2000 R² = 0.675 1000 0 0 0 1000 5000 2000 3000 4000 Station (mm) 5000 0 6000 3000 2000 R² = 0.771 1000 0 0 1000 1000 6000 2003 4000 TRMM (mm) TRMM (mm) 2005 4000 2000 3000 Station (mm) 4000 2000 3000 Station (mm) 4000 5000 2007 5000 4000 3000 2000 R² = 0.697 1000 0 5000 0 1000 2000 3000 4000 Station (mm) Figure 5.8 Comparison annual rainfall between TRMM rainfall and station rainfall in Thailand 170 5000 6000 Rainfall Parameters TRMM (Doy) 200 SRS 180 160 140 R² = 0.508 120 100 100 140 160 Station (Doy) 180 200 300 400 LRS ERS TRMM (Doy) TRMM (Doy) 120 350 300 R² = 0.844 250 250 200 R² = 0.668 150 100 200 200 250 300 Station (Doy) 350 100 400 150 200 Station (Doy) 250 Figure 5.9 Comparison rainfall parameters between TRMM rainfall and station rainfall in Thailand 171 300 5.5 Discussion and Conclusion 5.5.1 Spatio-temporal Variation of Vegetation Phenology and Rainfall Seasonality The results presented in this research demonstrate that satellite time-series data provide great opportunities to study not only regional vegetation variability but also internal climatic fluctuation. These data can identify the relationships between phenology and climate variability, as well as the drivers of phenological changes related to climate variations. The rainy season in this region tended to start early and end late resulting in a slightly longer length. However, the amount of rainfall has decreased from 2001 to 2010. Annual rainfall decreased across most of the western and northern regions of the Peninsula and increased across the eastern region of Thailand. The results reveal a significant positive correlation between SGS and SRS overall, although negative correlations can be found in areas used most extensively for agriculture. The forest areas and rainfed croplands show similar patterns, which shifted the growing season according to changes in rainfall, especially in Myanmar. The rate of phenological changes with respect to rainfall seasonality varied among different land cover types and ecosystems. In addition, the driver of phenological change due to rainfall variability and human management was also observed on a local scale. The local scale analysis indicates the hotspots and sensitive areas of vegetation conditions and climate variability. For example, the results showed that the northwestern part of the Peninsula was distinctly different from other regions in the Peninsula. A later start and earlier end of growing seasons with shorter length were found in this area. Moreover, the rainy season was later and shorter than it was in the rest of the region and the amount of rainfall in this area decreased over the 10 years study. This area is covered mainly by forests and has low human disturbance; as a result, phenological patterns 172 reflect the strong influence of climate in this location. This sensitivity of phenology to rainfall in some specific sites is important for environmental management. Additionally, the findings indicate that the major drivers of the phenological changes are climate and human management. The regional and local scale analyses indicate that rainfall is a dominant force in naturally vegetated areas and rainfed croplands, whereas human management is a key factor in heavily agricultural areas with irrigated systems. 5.5.2 The Influence of Climate Change on Phenological Patterns The general pattern and phenological changes exhibited relative correspondence with rainfall seasonality in SEA. The wet and dry years occur as a repeated cycle in this region. Extreme climate events in this region are associated with the El Niño Southern Oscillation Phenomenon (ENSO) (NIC, 2009). ENSO is related to the periodic shift in global climate resulting in a pattern of seasonal monsoons. Monsoon frequency and intensity enhanced by ENSO can create or intensify floods and droughts (NIC, 2009). According to Ninh & Kelly, 2000, El Niño and La Niña are used to describe the periodic warming and cooling of the tropical Pacific Ocean and the consequent disruption of the atmospheric circulation, bringing extreme weather and climate to many low-latitude areas. Both El Niño and La Niña events have severe impacts on the Indochina region, affecting patterns of temperature, rainfall and other weather variables such as the frequency of tropical storms. An El Niño event is associated with drought in most parts of the region, whereas La Niña events are more likely to bring excessive monsoon rain leading to flooding, especially in low lying areas. SEA experiences an increase of rainfall of up to 150% and frequent flash floods have been reported in SEA related to La Niña events (ADB, 2009). 173 When considering the EVI and phenological patterns in this research, it was apparent that the change of EVI and phenology were evident as an irregularly repeated cycle corresponding to the years of El Niño and La Niña. The EVI and rainfall showed decreasing patterns in 2002, 2004-2005, 2009-2010 (Figure 4.7 and Figure 5.3). This pattern was related to El Niño years with very low rainfall (TMD Thailand, 2010), resulting in lower EVI values than those in other years. La Niña years, which are wet with high amounts of rainfall, occurred in 2003 and 2007-2008 (TMD Thailand, 2010). EVI in these years illustrated larger increases, particularly in 2007-2008. In addition, phenological patterns also revealed the change and shift patterns during those extremely climatic years (Figure 4.16). Therefore the temporal EVI and rainfall, including phenological patterns showed the cycles of changes (e.g., decreasing, increasing, change, and shift patterns) over the ten year study period, especially with significant changes in the years of extreme climate events. Several reports revealed the impact of extreme climate events on vegetation condition in this region. These extreme years have caused floods during the rainy season and water shortages during the dry season. For example, the model of Naylor et al., as cited in NIC, 2009) predicted that El Niño events typically cause a delay of up to 2 months in the onset of the monsoon season in Indonesia, which delays the planting of rice crops. Lansigan et al., (as cited in NIC, 2009) investigated the potential impacts of climate during El Niño and non-El Niño years in the Philippines. They found a major change in the cycle of planting and harvest during El Niño years: the variability in weather patterns moved the sowing date for rice to as early as Julian day 137 (mid-May) or as late as Julian date 229 (mid-August) when it is normally near Julian day 173 (mid-June). Additionally, changes in precipitation patterns affected current cropping patterns, crop growing seasons, and the sowing period (ABD, 2009). 174 In addition to the growing season shift, the precipitation variability in El Niño years reduced crop yield. An El Niño year was typically marked by a shorter and more intense wet season, which decreased crop growing time and subjected crops to stress from excess water. This situation caused a crop loss of 52-81 percent, and the delayed onset of the rainy season by up to a month in ENSO years cut rice yields by up to 11 percent in parts of Indonesia (NIC, 2009). Boonpragob (2005) revealed that the economic losses in Thailand between 1989 and 2002 due to floods, storms, and droughts were mainly from the agriculture sector because crop yield losses amounted to more than $1.25 billion. In recent years, thousands of hectares devoted to rice production have been damaged by frequent flooding in the Red River Delta, Central Region, and Mekong Delta in Vietnam due to the impact of extreme weather events (ADB, 2009). The extreme climate events are particularly important for SEA because the region’s economics depends on agriculture (Ninh & Kelly, 2000). Agricultural outputs are severely affected if there is a serious deficit in water supply or excessive rainfall. Therefore, climate change has been and will continue to be a critical factor affecting the productivity of agriculture in this region (Lasco et al., 2011). In addition, the increase in incidence of observed climate extremes such as floods, droughts and tropical cyclones in this region have caused adjustments in farm management practices in response to climatic variation. The most commonly used adaptation techniques in the agriculture sector of SEA involve changes in crop variety, cropping patterns, cropping calendar, and improved farm management (Lasco et al., 2011). As a result, these adjustments lead to changes in phenology in this region; for example due the later start of rainy season, farmers have to delay the growing season. Therefore, climate variability will continue to cause significant impacts on SEA agriculture in coming decades. The frequency and severity of climate extremes due to climate change are expected to cause greater losses, 175 which will be exacerbated if extremes occur during vital stages of crop growth. Most importantly, these impacts are critical to food shortages in the region and pose a serious threat to future food security. Not only has climate variability had substantial impacts on agriculture, but it has also negatively affected naturally occurring vegetation. Changes in forestry and vegetation due to climate change will likely impact species biodiversity. The biologically diverse forests of SEA represented 5.1 percent of the total forest areas in the world in 2005, and they were a major source of global forest products, accounting for 50 percent of total forestry exports from Asia and the Pacific (NIC, 2009). SEA forests are vulnerable to climate change due to degradation and unsustainable practices, such as illegal logging and conversion of native forests to agricultural lands. The report in Thailand (Boonpragob and Santisirisomboon, as cited in NIC, 2009) indicated that increases in temperature and variations in precipitation in Thailand due to climate change will cause the expansion of tropical dry forest into subtropical moist forest in the northern part of the country, and cause replacement of subtropical forests with tropical forests in the southern part of the country. In addition, the report in Vietnam concluded that climate change is likely to cause significant alterations in the Vietnamese forestry industry, particularly to shifts in suitable growing regions (NIC, 2009). These reports suggest that climate change poses a threat to biodiversity in the region. Subsequently, climate change will cause the replacement of highquality forests with low-quality forests, which will likely to lead to significant biodiversity loss. Thailand and Vietnam are expected to suffer more than other countries from biodiversity loss due to the impact of future climate change (ADB 2009). EVI and phenological patterns including rainfall seasonality derived from satellite images provide an opportunity to study phenological information related to climate changes over 176 large spatial regions and enable environmental monitoring and assessment. As observed from the results of this research, phenological patterns and changes can identify relationships with climate change events, particularly in El Niño and La Niña years. Additionally, the results in a local scale analysis can identify the hot spots or sensitive areas of decreasing EVI and forest loss as well as significant changes of phenology in croplands. These results could guide management priorities for protection and conservation in forest areas and agricultural management in cropland areas. SEA is more vulnerable to the impacts of climate change with its rapidly-growing population and increasing dependence on natural resources and agriculture (ADB, 2009; Lasco et al., 2011). Consequently, the information about the relationship between phenology and climate patterns is critical for examining the effects of climate change on ecosystem dynamics in this region on regional and local scales. 5.5.3 The Influence of Human Management on Phenological Patterns In addition to the influence of climate on phenological changes, land use is also an important determinant in this region. Agricultural land use is one of the most extensive land cover types in this region. The spatio-temporal patterns of phenology in agricultural areas are distinct from that of natural vegetation. This characteristic can be seen at both regional and local scales. This study shows that the trends of phenological parameters reflect changes in land use practices in primary agricultural areas. However, the timing and production of crop phenology can be quite variable and dependent on crop types and human management. Socio-economics and improved farm management are the main causes of land use change in this region. The increase or decrease in crop prices would influence changing in crop types by farmers. In addition, changes in policy and planning could result changes in the extent and intensity of cultivation. These factors affect the timing and production of crop phenology to 177 change over time. Nevertheless, a mixture of climatic and anthropogenic effects is also possible. Extreme climate events such as floods and droughts have also resulted in the land use changes, for example, changes in cropping patterns (ADB, 2009; Lasco et al., 2011). More data on these changes are required in order to perform for a site by site analysis. The results presented in this research indicate that some locations are primarily influenced by human management, but climate also plays a role in some locations. However, the phenological patterns and changes that are influenced by human management are likely to become an increasingly important issue in this region. 178 CHAPTER 6 ENVIRONMENTAL IMPACTS WITH BIOGEOCHEMICAL MODEL: RICE FIELDS IN THAILAND 6.1 Introduction Chapters 4 and 5 examined spatio-temporal variations of vegetation phenology and rainfall seasonality, and responses of seasonal vegetation dynamics to climate variations in SEA. In general, phenological changes in SEA exhibited relative correspondence with rainfall seasonality and major climate events. However, human management with land use practices also altered land surface phenology in some areas. Therefore, the major drivers of the phenological changes were both climate and human management. To better understand the environmental consequences of climate change and management practices, it is necessary to assess the impact of environmental change, particularly in agricultural areas. This chapter demonstrates the integration of remote sensing mapping with a biogeochemical model (DNDC) to quantify GHG emissions (CO2, CH4, and N2O) from rice fields in Thailand. Lopburi province in Thailand was selected as a case study because this site is an important rice field site in the Chaopraya River Delta in Thailand. The new approaches described in chapter 3 that include multi-temporal remote sensing data and the new database system with a grid-based unit at 250 m resolution, were proposed to improve model performances and provide more accurate emissions predictions. This approach can simulate emissions at site level, which reflects the spatial diversity of crop growth environments, and can also advance the regional estimation of the DNDC model. The database was established from climate and soil data, in combination with additional farm management information to estimate emissions for single and double rice cropping systems. 179 Additionally, spatio-temporal patterns and changes of GHG emissions were investigated under different scenarios (actual data, different rates of fertilizer, climate variation, phenological changes, and climate variation and phenological changes) to identify the effects of the input parameters on the model regional emission (see more details in chapter 3). The emission rates were directly retrieved cell-by-cell from the model, while total yearly emissions were obtained by multiplying the emissions by the area of the pixels containing rice fields (multiply by the size of the grid cell if that cell has rice fields) to produce total emissions per cell for scenario1 (actual input data). The regional emissions of the province were derived from a sum of all cells in the province for simulated year. To indicate the net effect on global warming, Global Warming Potential (GWP) was calculated for each grid cell. Furthermore, comparisons of the DNDC results to the IPCC approach and Thailand research is performed in this chapter. The influence of human management, climate variation, and physical geography on the change of GHG emissions is also discussed. 6.2 Spatio-temporal Patterns of GHG Emissions under Different Scenarios This research conducted five scenarios for a selected site (Table 3.5) to identify the effect of management practices (fertilization, flooding period), climate, and phenology on modelproduced regional emissions. The results of five scenarios are explained in the following: 6.2.1 Scenario 1: Actual Input Data for 2002 and 2010 Tables 6.1-6.2 and Figures 6.1-6.3 illustrate the results from the regional simulations based on actual input data (climate, land use, phenology, farming management) in Lopburi Province in 2002 and 2010. The overall emissions results at this site indicate that all emissions increased in 2010. The modeled average of CH4 for the simulated years of 2002 and 2010 were 206 and 255 kg CH4-C/ha, respectively. The relative deviation of CH4 emission rate 180 in 2010 was higher than 2002 at 24%. The average of CO2 was 725 kg C/ha in 2002 and 770 kg C/ha in 2010 with a relative deviation of 6%. The average of N2O was 14 kg N/ha in 2002 and 17 kg N/ha in 2010 with a relative deviation of 16%. Yearly total emissions and GWP present a similar pattern, where emissions were higher in 2010 than in 2002. The sum of the warming forces of all the three GHG emissions based on the concept of GWP shows that GWP values increased 18% from 2002 to 2010. The total emissions of CH4 in this province were 43,707 tons CH4-C in 2002 and 61,729 tons CH4-C in 2010. The total emissions of CO2 were 154,152 tons C and 186,563 tons C in 2002 and 2010, respectively. The total emissions of N2O in 2002 were 3,055 tons N and 4,042 tons N in 2010. The simulated spatial distribution of GHG emissions from the Lopburi rice fields and their GWP under continuous flooding condition is shown in Figures 6.1-6.2 The higher emission rates of CH4 were found in the western region of Lopburi province (Amphoe Ban Mi and Tha Wung), which is dominated by a double cropping rice system (Figure 3.9). The lower CH4 emission rate occurred in the central area of Lopburi province (Amphoe Khok Samrong), which is dominated by a single cropping rice system. It can be noted that CH4 showed sensitivity to soil properties. The patterns of CH4 reflected spatial variation in the soil properties (soil texture) when comparing the result with the soil maps (Figure 6.4). Low CH4 emissions were located in soils with a high clay fraction (higher than 50%). Conversely, high CH4 emission rates were found in low-clay soils, even in single cropping rice systems (yellow and orange colors in Figure 6.1). 181 However, CO2 and N2O displayed different spatial patterns from CH4. Emission rates of CO2 and N2O showed high values in single cropping systems and low values in double cropping systems. This is due to several influencing factors. Low-clay soil in areas with a single rice cropping system produces more CO2 and N2O emissions. In addition, a single cropping system has less crop biomass production than a double cropping system leading to a greater loss of soil organic carbon that affects higher CO2 emissions. N2O emissions were lower during the flooding period; as a result, shorter flooding period lengths for a single rice cropping system could produce more N2O emissions. High emission rates were found in Amphoe Sa Bot for all GHG emissions mainly due to the soil properties. When comparing the results of emission rates in 2002 with 2010, CH4 clearly showed an increasing pattern in a double cropping rice system. CO2 also showed increasing emissions in a double cropping rice system; however, decreasing emissions were found in a single cropping system. N2O emissions increased in most areas of the study site, although the area with the highest N2O emissions (orange color in Figure 6.1) in 2002 showed a slight decrease in emissions in 2010. Table 6.1 Emission rate in scenario 1 Parameter Emission Rate 2002 2010 % Change CH4 (kg CH4-C/ha) 205.69 254.62 23.79 CO2 (kg C/ha) 725.44 769.53 6.08 N2O (kg N/ha) 14.38 16.67 15.96 16,249.03 19,116.33 17.65 GWP (kg CO2-equiv/ha) 182 Table 6.2 Total emission in scenario 1 Parameter CH4 (Ton CH4-C) Total Emission 2002 2010 % change 43,706.85 61,728.50 41.23 CO2 (Ton C/ha) 154,151.81 186,562.70 21.03 N2O (Ton N/ha) 3,055.19 4,042.05 32.30 3,452,817.86 4,634,515.16 34.22 GWP (Ton CO2-equiv/ha) 6.2.2 Scenario 2: Fertilizer Effect This scenario used different fertilizer rates for 2002 and 2010 (with a higher rate of fertilizer in 2010) in order to test the effect of fertilization on GHG emissions. Tables 6.3 and 6.4 present the overall scenario 2 results and Figures 6.3 and 6.5 show spatial and temporal patterns of GHG emission rates and GWP. The results from the regional simulations based on different fertilizer rates in scenario 2 indicate that CH4 and N2O emissions in 2010 were higher than 2002; however, CO2 presented the opposite pattern with lower emissions in 2010. The mean emission rate of CH4 increased by 23% from 2002 to 2010 (Figure 6.3). A decrease of 5% was found for CO2. The mean emission rate for N2O showed only a 1% increase from 2002 to 2010. These results suggest that a change in fertilizer rate can significantly affect CH4 emissions. CH4 demonstrated a higher rate of increase, which was similar to scenario 1 (Figure 6.4). On the other hand, N2O emissions in scenarios 1 and 2 were quite different asN2O emission rates showed an increase of 16% in scenario 1. GWP showed a slight increase at 9.7%. 183 Figure 6.1 GHG emission rate in scenario 1 184 Figure 6.1 (cont’d) 185 Figure 6.2 Total GHG emission in scenario 1 186 Figure 6.2 (cont’d) 187 CH4 1000 CO2 300 900 kg C /ha kg CH4-C /ha 250 200 150 S1 100 S2 700 S1 600 S2 500 50 0 400 2002 2010 2002 2010 5000 20 kg CO2 equi.v /ha N20 15 kg N/ha 800 10 S1 5 S2 0 GWP 4000 3000 2000 S1 1000 S2 0 2002 2010 2002 2010 Figure 6.3 Comparison GHG emission rate between scenario 1 (S1) and scenario2 (S2) 188 Figure 6.4 Soil texture and % clay in Lopburi province 189 Figure 6.5 clearly illustrates the spatial variations of emissions across the province at grid-cell scale. CH4 emissions increased in 2010, particularly in double rice cropping systems. In contrast, a decrease in CO2 emissions was found in a double rice cropping system. Increases in N2O emissions were sparsely distributed in a double rice cropping system. Therefore, the results of this scenario indicate a strong effect of fertilization on CH4 emissions. In addition, increases in fertilizer rate showed a significant impact on a double rice cropping system, which exhibited large increasing values of emission, due to the increasing frequencies of fertilizer application in this cropping system. Table 6.3 Emission rate in Scenario 2 Parameter Emission Rate 2002 2010 % Change CH4 (kg CH4-C/ha) 205.69 253.90 23.44 CO2 (kg C/ha) 725.44 690.60 -4.80 N2O (kg N/ha) 14.38 14.57 1.35 16,249.03 17,819.65 9.67 GWP (kg CO2-equiv/ha) Table 6.4 Total Emission in Scenario 2 Parameter CH4 (Ton CH4-C) Total Emission 2002 2010 % change 43,706.85 53,952.91 23.44 CO2 (Ton C/ha) 154,151.81 146,749.17 -4.80 N2O (Ton N/ha) 3,055.19 3,096.51 1.35 GWP (Ton CO2-equiv/ha) 3,452,817.86 3,786,563.94 9.67 190 6.2.3 Scenario 3: Climate Effect The objective of this scenario is to investigate the effects of climate variations on GHG emissions by using actual climate data (temperature and rainfall) for nine years, while keeping all other input parameters constant (used input data from 2002 for all years from 2002-2010). Figure 6.6 shows a graph of emission changes. The spatial patterns of GHG emissions from 2002 to 2010 are included in Appendix D. The results illustrate little variability in CH4 emission during the nine-year period, except in 2010, which showed a significant decrease in emissions. CO2 also demonstrated a stable pattern but exhibited a slight decrease in 2007. On the other hand, N2O emission patterns varied across the nine-year period. When comparing with the climate input data (Figure 6.6), it can be noted that N2O followed the pattern of rainfall. There is also an indication that N2O increased with an increase of rainfall. GWP patterns were similar to N2O. In addition, the decrease of CH4 and CO2 in 2010 and 2007 were largely related to extreme climate events in this study area. When examining the climate data at this site, the temperature was generally higher and the rainy season arrived later with a shorter length in 2010 compared to other years; thus, 2010 was the drought year. In addition, the length of rainy season was longer with very low rainfall in 2007. This extreme climate with changes in temperature and rainfall could affect the change in CO2 and CH4 emissions. The results of this scenario suggest that rainfall could be a major factor for N2O emission and extreme climate events have an impact on CH4 and CO2 emissions in this region. 191 Figure 6.5 GHG emission rate in scenario 2 192 Figure 6.5 (cont’d) 193 6.2.4 Scenario 4: Phenology Effect This scenario examines the effects of phenological changes on GHG emissions. The actual phenology (the start and end dates of the growing season) of nine years were utilized in this scenario, whereas all other inputs were held constant by using data from 2002 for all years. The temporal CH4 emission varied greatly in this scenario (Figure 6.7). CH4 emissions showed decreasing patterns from 2002 to 2005 but exhibited a large increasing trend from 2009 to 2010. The temporal CO2 emissions illustrated increasing trends with a remarkable increase in 2005. N2O emissions demonstrated an increasing trend with small variations. Areas with the highest emissions for CH4 and CO2 did not change much during the nine-year period. Spatial patterns of GHG emissions for this scenario are provided in Appendix D. The modeled results indicate that phenological changes have significant impacts on emissions, particularly for CH4 emissions. Figure 6.7 compares CH4 emissions with the length of growing season; these two graphs show very similar patterns, as the length of growing season increased, CH4 increased. This is due to the fact that the flood duration for rice paddies is related to the length of growing period, especially in a double rice cropping system. A longer growing season length needs a longer flood duration (no mid-season drainage is included in the simulation of this research). The flooded rice paddy provided a favorable environment for methanogenesis, resulting in increased CH4 emissions. Therefore, change in the number of growing seasons from a single to double rice cropping system affects CH4 emissions due to increases in the frequency and duration of the flooding period (e.g., twice in a year) in the paddy soils. 194 In the case of CO2 emissions, change in emissions was in relative agreement with vegetation production. Small integral and EVI, which shows the vegetation production or biomass, demonstrated decreasing trends from 2003 to 2005, but CO2 showed an increasing trend during these three years. This result is largely related to the loss of soil organic carbon. The decrease in crop biomass production or vegetation growth leads to the changes in soil decomposition rates and crop residue incorporation into the soil. This condition enables the change in soil from a sink to a source of atmospheric CO2 resulting in the increase in CO2 emissions. N2O showed very little changes in this scenario and its pattern was similar to GWP. 6.2.5 Scenario 5: Climate and Phenology Effect The actual climate and phenological data for the nine-year period were utilized in this scenario and the remaining input parameters used were from the 2002 dataset. The results show a significant impact of climate and phenology on annual emission rates (Figure 6.8). Spatial patterns of GHG emissions of this scenario are shown in Appendix D. The trends of all emissions were similar to the trends of scenario 4. In the case of CH4 emissions, there was a decreasing trend in the first four years (2002-2005) and an increasing trend from 2009 to 2010. CO2 and N2O emissions in this scenario have more variation when comparing to the trend of scenario 4. Both emissions illustrated rapid changes in 2005. These results suggest that climate combined with phenological changes have significant influences on CH4 emission and also affect the variation in CO2 and N2O, particularly during the years of extreme climate events in 2004 and 2005. 195 Emission Rate 35 CH4 230 210 33 190 32 2010 2009 2008 2007 25 CO2 800 600 Mean Min Temperature 24 23 °C 400 22 200 21 0 20 20 1600 15 1200 MM N2O 10 800 5 400 0 0 200 20000 GWP 18000 Lenght of Rainy Season 160 Days 16000 14000 120 80 2010 2009 2008 2007 2006 2005 2004 2002 2010 2009 2008 2007 2006 0 2005 10000 2004 40 2003 12000 2002 Kg CO2 equiv./ha Mean Rainfall 2003 Kg C/ha 1000 2006 2005 2004 30 2003 31 2002 170 150 Kg N/ha Mean Max Temperature 34 °C Kg-CH4-C/ha 250 Climate Input Data Figure 6.6 Comparison between climate input data and GHG emission rate in scenario 3 196 Emission Rate 350 CH4 Kg-CH4-C/ha Kg-CH4-C/ha 250 Phenology Input Data 230 210 190 170 Kg C/ha 1000 200 0.50 0.45 0.40 0.35 0.30 0.25 0.20 CO2 800 600 400 200 0 20 Kg N/ha 250 150 150 50000 15 EVI Small Integral 40000 N2O 10 30000 5 20000 10000 0 20000 GWP 18000 16000 14000 12000 2010 2009 2008 2007 2006 2005 2004 2003 10000 2002 Kg CO2 equiv./ha Length of Growing Season 300 Figure 6.7 Comparison between phenology input data and GHG emission rate in scenario 4 197 When comparing the results of this scenario with scenarios 3 and 4, it can be noticed that phenology is the primary factor effecting the changes in CH4 emissions. The graph of temporal CH4 emissions shows large variation with phenological changes in scenario 4. CO2 emissions trends were stable in scenario 3, 4, and 5 but showed rapid change in the extreme climate years. However, N2O emissions were strongly related to climate variation, especially rainfall changes. 1000 CH4 230 210 190 600 400 200 0 20 N2O 15 10 5 Kg CO2 equiv./ha 170 150 Kg N/ha CO2 800 Kg C/ha Kg-CH4-C/ha 250 0 20000 GWP 18000 16000 14000 12000 10000 Figure 6.8 GHG emission rate in scenario 5 6.2.6 Global Warming Potential of Lopburi Province Different GHG emissions can be directly compared by converting the fluxes of the non-CO2 into CO2 equivalents via their radiative forcings (IPCC, 2007). The radiative forcing of each gas was determined using its heat trapping capacity as well as lifetime in the 198 atmosphere. This research applied the GWP standard of the IPCC approach from 2007. Although the results indicate that the contribution of GWP over all scenarios was positive and varied from year to year, GWP increased in scenarios 1, 2, and 4, decreased in scenario 3, and was fairly stable in scenario 5 (Figures 6.3, 6.6-6.8). GWP showed fluctuant patterns with similarity to N2O in all scenarios. This implies that GWP is highly dependent on N2O according to the higher weight of the warming forces of N2O, compared to other emissions, in the GWP equation. The spatial patterns of GWP indicate that areas with high GWP values also have a high rate of N2O. However, CH4 and CO2 also shares in GWP especially in the areas with high emission rates of these two gases. These results suggest that reducing N2O could effectively mitigate the net effect of rice fields at this site on global warming. Additionally, the positive values of GWP imply that, in general, paddy rice fields are a source of atmospheric GHG at this study site. 6.3 Comparison to IPCC Estimation and Research in Thailand A baseline emission factor from IPCC guidelines (IPCC, 2006) is usually applied for estimating the regional or global CH4 emissions from rice paddies, but none for CO2 and N2O. With the IPCC emission factor of 130 CH4-C/ha (Table 6.5) under 100-day continuously flooded conditions without organic amendments for rice paddies, the total CH4 emissions of Lopburi province were 27,624 tons CH4-C in 2002 and 31,517 tons CH4-C in 2010. The modeled estimation of regional CH4 emissions in this study site from the actual data (scenario 1) was 43,707 and 61,729 tons CH4-C in 2002 and 2010, respectively, and the mean emission rate was 199 206 and 255 CH4-C/ha in 2002 and 2010, respectively. These results indicate that CH4 emissions from the DNDC model were greater than that of the IPCC approach. The relative deviation was found to be -58% and -96% in 2002 and 2010, respectively. However, if the continuous flooding period changed from 100 days to 200 days due to a double rice cropping system, which is the major rice system in this site, the results demonstrate good agreement between modeled emissions and the IPCC approach. The emissions from the IPCC method (260 CH4-C/ha) were larger than the model estimates in 2002 (206 CH4-C/ha) but were very close to model estimation in 2010 (256 CH4-C/ha). The relative deviation was only 26% and 2% in 2002 and 2010, respectively. The CH4 emission factors and scaling factors from previous Thailand research (Gale et al., 2005) were also used to compare with the regional CH4 emissions of this research. Gale et al. (2005) reported both the standard scaling factors for Thailand (232 CH4-C/ha) and the CH4 emissions rate for each province of Thailand in 1998 (Table 6.6). The standard estimate rate of the Thailand research showed close agreement with the modeled estimation. A relative deviation between the standard and the modeled estimation was only 12% and 10% in 2002 and 2010, respectively. In addition, estimated CH4 emission rates for each province in Thailand based on minimum, median, and maximum scenarios in 1998 (Gale et al., 2005) were also compared to modeled emission rates (Table 6.7). The results demonstrated that the modeled emission rate was in between the median and maximum scenarios (150-331 CH4-C/ha). 200 Therefore, the results indicate close agreement between the standard emissions rate and modeled emissions rate with the IPCC approach under 200-day continuous flooding and scaling factors from the Thailand research. The discrepancies between simulated emissions, emission factors and standard emission rates in the IPCC approach and the Thailand research may be due to the differences in climate and cropping practices in each simulated year. The DNDC model has a number of advantages over the standard scaling factor method. The DNDC model is a complex, process-based model that applies various ecological drivers and soil environmental variables for regional GHG emissions. The greater detail of the DNDC calculations allows the effect of different soils, crops, climates, and farming practices to be reflected in the emission estimate. Climate and farming management can vary significantly in different years. The DNDC model is able to vary the input parameters (climate, soil, vegetation, and management factors) for different years and different grid cells. This technique has the potential to greatly improve the accuracy and precision of GHG estimates and reduce uncertainties in estimated parameters concerning rice field management. In addition, the DNDC model allows users to identify and assess the specific effects of any single model driver or model input parameter on GHG emissions. This approach may be applied from local to regional scales. Table 6.5 Comparison between modeled estimation and IPCC approach DNDC IPCC 2006 (100-day continuously flooded) Year IPCC 2006 (200-day continuously flooded) Kg CH4Kg CH4Kg CH4Ton C Ton C RD (%) Ton C RD (%) C/ha C/ha C/ha 2002 205.7 43,706.8 130.0 27,624.2 -58.2 260.0 55,248.4 26.4 2010 254.6 61,728.5 130.0 31,516.9 -95.9 260.0 63,033.8 2.1 RD = relative deviation Source: IPCC. (2006). Guidelines for National Greenhouse Gas Inventories, Volume 4: Agriculture, Forestry and other Land Use. (Emission Factor: 1.30 kg CH4/ha/day) 201 Table 6.6 Comparison between modeled estimation and Thailand standard scaling factors Year DNDC Thailand research1 (2005) RD (%) Kg CH4-C/ha 232.4 11.5 232.4 -9.6 Kg CH4-C/ha 2002 205.7 2010 254.6 RD = relative deviation Source: Gale et al. (2005). Development of a database for estimating methane emissions from rice fields in Thailand (Emission Factor: 232.36 kg CH4/ha) Table 6.7 Comparison between modeled estimation and standard emission rate in Lopburi province Year DNDC Kg CH4-C/ha 2002 205.7 Thailand research2 (1998) Min RD Median RD Max RD (Kg CH4C/ha) (%) (Kg CH4C/ha) (%) (Kg CH4C/ha) (%) 55.1 -273.6 150.1 -37.1 330.8 37.8 2010 254.6 -362.5 -69.7 23.0 RD = relative deviation Source: Gale et al. (2005). Development of a database for estimating methane emissions from rice fields in Thailand (Estimated CH4 emission rate for each province in Thailand based on min, median, and max scenarios) 6.4 Discussion and Conclusions This research demonstrates the potential integration of remote sensing mapping with a biogeochemical model to quantify spatio-temporal patterns of GHG emissions from rice fields in Thailand. Spatial and temporal dynamics of GHG emissions were characterized under different scenarios to investigate the effects of various factors on GHG emissions from rice fields. These factors included climate and human management (fertilization and phenology). The results of five scenarios estimating GHG emissions from rice fields in Thailand demonstrate the influence of human management, climate variation, and physical geography on the change in GHG emissions. The results of GHG estimations indicate that phenology is the main factor affecting 202 the changes in CH4 emissions. The change of CO2 emissions was relatively smaller than CH4 and N2O in all scenarios but showed rapid changes in extreme climate years. N2O emissions were strongly related to climate variation, especially rainfall changes. Additionally, high CH4 emissions showed a correspondence with light texture soil. An increase in rates and applications of fertilizer produced the high CH4 and N2O emission. A longer flooding period resulted in increased CH4 emissions. 6.4.1 Effect of Soil Properties on GHG Emissions Although GHG emissions in rice fields are influenced by a variety of factors, soil properties are some of the most significant factors in changing GHG emissions. The results indicate that the chemical and physical soil factors played important roles in CH4 and N2O emissions. The soil chemical properties, soil organic carbon, and soil texture, affected CH4 emissions. Emissions increased with as clay content decreased and soil organic carbon increased. This finding is in agreement with numerous observations (Li et al., 2004, Babu, et al., 2006, Horwath, 2011, Smakgahn, 2003). High organic content and light soil texture contain abundant carbon sources for methanogenesis which leads to high emission rates. In contrast, heavily textured soils (high percent of clay) lack the nutrition for methanogenesis and emissions are entrapped before they are released to the atmosphere (Smakgahn, 2003). Additionally, heavier textured soils emitted less CH4 than lighter soils under the same management conditions due to the clay absorption which limits DOC availability to the soil microbes (Li et al., 2004). Horwath (2011) also mentioned that the higher buffering capacity of clay soils can be directly related to 203 lower CH4 emissions, which is due to the increased content of alternate electron acceptors resulting in a prolonged onset to reduced conditions. In addition, soils with a rapid return of pH to neutral upon flooding decreased in soil redox potential (Soil Eh) more rapidly and hastened the onset of CH4 production (Horwath, 2011). The same pattern is also found for CO2 and N2O emissions; heavier soils emitted less CO2 and N2O than lighter soils. In addition to soil chemical properties, soil physical properties also influence CH4 production and emission. Soil Eh is related to soil moisture level, with high soil moisture resulting in a low soil Eh. A decrease in soil Eh under continuous flood conditions leads to an increase of CH4 production rates. According to Smakgahn (2003), a low soil Eh leads to low oxygen levels in the soil, affecting the growth of methanogenic populations and CH4 emissions, whereas high soil Eh from dry conditions is not suitable for methanogenesis activity. The range of soil Eh for CH4 emission is −150 to − 210 mV. On the contrary, soil Eh for significant N2O emissions ranges between +120 to +250 mV (Smakgahn, 2003). When rice fields are drained, the result is low soil moisture with a high soil Eh. These conditions decrease CH4 emissions, but intensify N2O emissions, particularly when N fertilizer is applied as top dressing. 6.4.2 Effect of Management Practices on GHG Emissions In addition to soil properties, field management practices exhibited significant impacts on CH4 emissions. Although phenology may not be directly influencing the change in emissions, it is related to management practices resulting in the change in CH4 emissions. Water management strongly affects CH4 emissions. Flooding duration of rice fields depends on the 204 length of the growing period. Flooded rice with high cropping intensity (more than one season in a year) is extensively practiced in SEA due to their high rice demand. The longer the flooding period during a given year, the higher CH4 emissions are found in these rice flooded fields (Li et al., 2004, 2005; Pathak et al., 2005). For example, a double rice cropping system is usually flooded from the first day of planting until 15 days before harvesting. CH 4 is produced through anaerobic decomposition of organic matter in biological systems. CH4 is emitted during flooding because anaerobic conditions are essential for methanogenic bacteria, which are a source of methane. Wetland rice during flooding is consequently favorable for both rice production and methane production. Moreover, low soil Eh from moist conditions is suitable for methanogenesis activity. Soil Eh is positive at the start of the growing season and gradually declines when the fields are flooded, and increases again when the fields are drained. Therefore, the change in length of growing season, resulting in the change of flooding duration, causes changes in CH4 emissions. Consequently, the cropping frequency and phenology of rice are significant factors affecting the change in CH4 emission. In contrast, flooded rice fields produce low N2O emissions, while the drainage of rice fields increases N2O emissions. Several researchers report high N2O emissions from the soil of rice fields that are drained during the growing season (Li et al., 2004, 2005; Pathak et al., 2005; Smakgahn, 2003). This is due to the higher rates of nitrification and denitrification than under continuously flooded conditions. Draining rice fields may create suitable O2 availability in the soil for N2O production as an intermediate product in either nitrification or denitrification, 205 while flooding may create strict anaerobic conditions and restrict N2O emissions (Li et al., 2004, 2005; Smakgahn, 2003). Therefore, water management is important for emission mitigation. Although mid-season drainage or pre-harvest drainage reduces CH4 and CO2 emissions, this practice significantly increases N2O emissions (Li et al., 2004, 2005; Gale et al., 2005). Similar research in Thailand, China, and India reported low N2O emission observed under flooded conditions (Li et al., 2005; Pathak et al., 2005; Smakgahn, 2003). Another important management practice for GHG emission is fertilization. The second scenario of this research demonstrates that the fertilization rate, particularly of urea, was positively correlated to CH4 emission. Furthermore, double rice cropping systems tended to produce more CH4 emissions than a single rice system because the applications and rates of fertilization increase to more than two applications per year. Smakgahn (2003) found that urea applied as a top-dressing fertilizer significantly increases CH4 emissions in Thai rice fields and higher nitrogen levels in fertilizer resulting in higher CH4 emission rates. The influence of N fertilizer on N2O emissions was also found in this research, with an increase in fertilizer rates resulting in increased N2O due to nitrification processes. Smakgahn (2003) also reported that urea can produce high N2O emissions. Fertilizer application directly affects rice growth (plant density, amounts of leaf). An increase in biomass and grain yield leads to increased emissions (Smakgahn, 2003). N fertilizer application increases crop biomass or yields, which could indirectly enhance the CH4 206 production (Zhang et al., 2009). Therefore, rice intensification, with a longer flooding period and increased fertilizer applications and rates, substantially enhance CH4 production. 6.4.3 Effect of Climate and Phenology on GHG Emissions The changes in climate variation (temperature and rainfall) and phenology (start date, end date, and length of growing season) have a large impact on GHG emissions from rice fields. Although CH4 and CO2 did not show major changes with climate variations, the extreme climate events (e.g., a drought year with high temperatures and low rainfall) had negative impacts on CH4 and CO2. Rainfall is an influential factor on N2O emission, as greater rainfall tends to produce higher N2O emissions. Wassmann and Dobermann (2006) also indicated that N2O is primarily emitted in pulses after fertilization and strong rainfalls. In addition, phenology change is associated with climate variation; phenology showed significant changes in 2005 and 2010, which were El Niño years. These phenological changes, such as a change in the small integral or length of growing season, have effects on crop production resulting in changes of CO2 and N2O emissions. Most importantly, phenological change (e.g., the length of growing season) due to the intention to increase yield or climate change (e.g., drought or flood years) demonstrated a significant influence on GHG emissions. In particular, for CH4, change in the length of the growing season leads to change in flooding duration resulting in the modification of CH4 emissions. It is interesting to note that the trend of GWP shows significant influences of N2O emissions on global warming. This is because the weight of the warming forces for N2O is 207 higher than the weight for CH4 and CO2. Additionally, all scenarios have positive GWP values; this suggests that paddy rice fields are a source of atmospheric GHG. As a result, to mitigate the net effect of rice fields on global warming, N2O emission should also be considered as an important factor. 6.4.4 Challenges and Limitations Challenges and limitations should be considered when applying the DNDC model. The first challenge is the quality of the input parameters. Several input datasets are required for DNDC model (i.e., climate, soil, land use, and farming practices) but some data are unavailable or there is a lack of information from field observations or data records; therefore, default values in the DNDC model or estimated values are applied in this research to estimate emissions, such as crop residual. For example, one important input parameter is rice variety. Although this research applied actual yields for 2002 and 2010 in scenario1, only one type of rice was applied in the model due to lack of information. Rice varieties and rice cultivars vary in their phenology and physiological features, and the DNDC model can differentiate between two or more rice varieties and cultivars by using several specific input parameters, such as biomass yield, grain:shoot:root biomass ratio, crop growth, and crop phenology. Additionally, to obtain inter-annual estimation from the model, yearly land use and rice ecosystem maps for different locations are also important for the DNDC to reduce uncertainty in emission estimates. Sufficient and high quality input data are important to obtain high estimation accuracy when using the DNDC model. The second challenge is that the effectiveness of the model depends on the study site. Continuous modification and calibration can enhance the performance of the DNDC model for estimating GHG emissions. Therefore, field data are necessary for both modification and 208 evaluation of model performances. Validating the model with the standard emission factor is unable to reflect the change of input parameters (e.g., climate, soil, management practices). As a result, discrepancies exist between the simulated results and the standard emission. Therefore, field observations are required to validate the model and improve the efficiency of model. This research generated phenology for both single and double rice cropping systems but it is possible that some pixels have no phenological data due to missing values or atmospheric effects. This condition might affect emission estimates, particularly for CH4 emission. This is because CH4 emission estimates requires start and end dates of the growing season to calculate flooding period. Water management is another important factor for CH4 and N2O emissions. This research applied only continuous flooding management. Different water management practices should be applied to obtain accurate estimations and to identify the effects of water management under different conditions. In addition, the change in GHG emissions reported in this research should be considered as the results of a specific site, which may not be representative of all climate/soil/management regimes across the world’s rice fields. However, the general trends presented in this research could be applicable to other locations, particularly in SEA. 6.4.5 Potential of DNDC Model The DNDC model integrates influences of various ecological drivers (climate, soil properties, vegetation, and human activities) and soil environmental variables for regional GHG emissions. In addition, the DNDC model can simulate daily emissions and soil Eh patterns throughout the rice-growing season and can define wet-dry conditions in flooded rice fields. 209 Additionally, the model provides graphic data for users to observe the patterns of emissions, and it is possible to estimate emissions on a large scale and generate mapped results. This data is useful for predicting the production potential and emissions from rice fields. The method introduced in this research, the DNDC in site mode with grid-based units to provide output on a regional scale, could play an important role in linking management changes to biogeochemical cycles based on spatially differentiated information. This method has the ability to preserve the advantages of site-based modeling while also meeting the demand for large-scale estimation. Model efficiency is enhanced by this technique. Therefore, the integration of remote sensing mapping with the DNDC model to quantify emissions can increase the accuracy of GHG emissions on a regional scale. Additionally, satellite derived phenology is a type of alternative information that can be utilized in the DNDC model because in-situ phenological information (e.g., start and end dates of the growing season) is difficult to obtain in some areas, and phenology from remote sensing could reduce the cost of ground survey. This method would greatly improve estimations compared to the IPCC method based on the baseline emission factor. Although the results of GHG estimations from the DNDC model need more validation data, the DNDC model can demonstrate patterns and trends of emissions and it is useful for monitoring emission changes and developing mitigation strategies. In conclusion, climate variation has both direct and indirect impacts on GHG emissions. Climate variation alters plant-soil systems resulting in GHG emission changes. In addition, climate change indirectly affects GHG emissions by influencing changes in phenology and management practices. Moreover, a phenological change is strongly related to management practices that lead to a change in GHG emissions. The results of this research indicate that GHG emission not only vary spatially but also change significantly over time. This research 210 demonstrates that extreme climate events (flood and drought) have a significant impact on CH4 and CO2 emissions. In addition, the most sensitive factor for N2O emissions from rice fields is the change in rainfall. The effects of climate variation also cause changes in phenology, which impacts emissions. For example, a change in the small integral (which represents seasonal vegetation growth) and the length of the growing season because of extreme climate events results in a change of emissions due to a crop production change. Human management in farming practices also produces huge environmental impacts. Two major management practices—flooding period and fertilization—affect emissions from paddy soils by altering the chemical and biophysical properties of the plant-soil system. Continuous flooding practices increase CH4 but lead to lower N2O emissions. An increase in the rates and applications of fertilizer enhances CH4 and N2O emissions. Additionally, soil properties influence GHG emissions, particularly for CH4. A light soil texture produces higher emissions than a heavy soil texture. Low soil Eh tends to increase CH4 emissions, but decreases N2O emissions. The results of this research suggest that practical mitigation options should be carefully regulated to balance the emission types more efficiently, as well as to maintain or improve grain yields. The results imply that mitigation of CH4 emissions could be emphasized for rice fields because its emission is much higher than CO2 and N2O. Soil Eh should be controlled by water management in the appropriate range (−100 mV and +200 mV) to reduce both CH4 and N2O emissions (Majumdar, 2003). However, it may be difficult to apply this management in areas with insufficient water supplies and during the rainy season. A reduction of N fertilizers can 211 reduce emissions but would also reduce total rice production. However, ammonium sulfate could reduce emissions by 10-67% (Gale et al., 200). The integration of remote sensing and environmental models can help in identifying strategies to reduce these GHG emissions while maintaining rice production, which is an important source of food in the world. 212 CHAPTER 7 THE INTERACTIVE PHENOLOGICAL ATLAS FOR SEA 7.1 Introduction Previous chapters illustrated spatio-temporal patterns and changes of phenology and rainfall seasonality as well as environmental assessment. This information is very useful for research of environmental and ecosystem dynamics. Due to the vast amount of information resulting from phenological and environmental analysis, a visualization system has become indispensable to represent these spatio-temporal processes and to simulate the dynamics of environmental systems in order to provide the understanding and knowledge of how these changes are linked to environmental consequences. This chapter presents the prototype of the Interactive Phenological Atlas (IPA) for SEA. This visualization system applies the concepts of visual exploration and Internet atlas to display thematic map, bivariate map, and map animation. IPA uses web mapping application to deliver phenological information and related environmental variables in SEA on Internet. The outputs from previous chapters were selected to be displayed in IPA. EVI, phenological parameters, rainfall seasonality, as well as annual data quality assessment maps were processed to represent spatial and temporal patterns of environmental changes in this region. Additionally, spatial distributions and changes of GHG emissions on a local scale analysis were also included in this system. The objective of this system is data exploration, which allows users to derive meaningful information from interactive processes. The target users are researchers in environmental and remote sensing fields. The application allows users to display and explore the information interactively. The basic functions in this system enable users to navigate the maps in both spatial 213 and temporal extents; moreover, dynamic graphs are performed over the advanced functions. Additionally, map animation provides the functions for users to explore, monitor, and compare patterns and trends of vegetation phenology in SEA. IPA provides opportunity for researchers to explore the spatial distributions, changes, and the relationships between phenology and climate variability, as well as GHG emissions in SEA. This information is important for the understanding of ecosystem dynamics and environmental consequences in this region. In this chapter, the major components in the design and implementation of IPA are described and the capabilities of IPA, including accessibility, data, navigation, and interactive functionalities, are explained. The chapter concludes with the achievement, challenges, and limitations of IPA. 7.2 Design and Implementation of IPA There are five major components in design and implementation of IPA which are explained in the following: 7.2.1 Graphical User Interface (GUI) Design The GUI for the main web page is the entrance to the IPA system that users can access with a web browser (Figure 7.1). This web page is designed to facilitate data exploration as well as to communicate information with simple accessibility to easy understanding of the web content without special experience or knowledge. There are six main panels for designing the GUI of IPA as shown in Figure 7.1. The details of these panels are described as follows: a) Map content: the map content covers the major area of the web page in order to enrich the information to users. The map elements in this section are the map, legend, navigation tools, based map overlay button, overview map, scale bar, and layer transparency adjustment. In 214 addition to the map elements, the dynamic time-series graph displaying the pixel-based information is also enclosed in this panel. This graph will be displayed as a popup window, when users chick on the specific location on the map. b) Layer content: there are two major map layer categories consisting of the general map layers (e.g., boundary, field points) and the thematic map layers (e.g., EVI, phenological parameters, rainfall seasonality, etc.). The checkbox is used to select any layers of the general map layers and the drop-down list allows user to choose any categories of the thematic map layers. c) Toolbar: this panel contains several tools for specific functions consisting of on/off legend, map animation display, and help functions. Users can select any themes of map animations from the drop-down list. Additionally, the latitude/longitude coordinate is also provided in this panel. 215 Figure 7.1 IPA interface and components 216 Figure 7.1 (cont’d) 217 d) Static graph display: there are two categories of graphs to be included in this section in order to show the temporal data for each environmental parameter. The first category is the biophysical parameters consisting of mean annual EVI and mean annual rainfall. Another category is the greenhouse gas emission including CH4, CO2, N2O, and GWP. e) Dynamic graph display: this small panel is overlaid on the map content to display the pixel-based time-series graph. The pop up window is displayed when user clicks on a specific location on the thematic map. f) Map animation: this panel displays on top of the map content panel when users click on the display animation button in the toolbar (Figure 7.1). The GUI of map animation provides interactive animation control and temporal legend. The start and stop time as well as duration of animation are set on the system. As a result, users can control animation by selecting play, pause, and stop buttons, and by selecting parameters on the top tool bar of the main web page. These tools allow users to explore, monitor, and compare the patterns and changes of phenology. 7.2.2 Map Design There are three categories of maps in this system—thematic map, bivariate map, and map animation. For thematic maps, the design is divided into two types for qualitative and quantitative data. The qualitative data, such as the start and end dates of the growing season, were designed by considering change of hue. The quantitative data applied the difference in lightness and saturation. The color scheme was designed according to data types. Sequential schemes were used for quantitative data such as the length of the growing season, whereas diverging scheme was designed for trend maps, which have negative and positive values. Sequential data classes are logically arranged from high to low values by sequential lightness 218 steps and diverging scheme emphasizes the values that progress outward from a critical midpoint of the data range. In addition to color scheme, data classification is also important. In this system, appropriate methods were selected for each map. Five classes were applied to most of the maps to reduce cognitive overload. However, some maps have more than five classes, such as the date of growing season that is related to the month of a year. Furthermore, the range of each class was assigned depending on the range of the raw data; for example, the range of trend map is the week corresponding to MODIS date. The bivariate map is required the specific design that are different from other maps. The bivariate map displays two attributes that allows users to examine and compare the relation between these two attributes. The important design is the selection of appropriate attributes and the use of a legend. This research created three bivariate maps (the start and end dates, length and amplitude, and length and large integral) by selecting two related phenological parameters in order to present the relationship of the two and to provide meaningful conclusions (Figure 4.16-4.18). Additionally, to create effective bivariate mapping, mean deviation technique was generated. Each map was produced by calculating the difference of phenological values of each year with the 10-year mean values based on pixel by pixel. The start and end dates represent the change in the timing of the growing season (early, normal, late). The amplitude, the length, and the large integral demonstrate the change of the growing season related to the annual biomass and net primary production (increase, normal, decrease). This information is able to indicate the sensitive areas of changes in vegetation condition as well as in land use. In addition to data selection, the legend design is another concern of this map. This system applied diverging/diverging schemes for bivariate maps because these maps use mean deviation technique to present normal, positive, and negative ranges. The technique was applied from 219 Brewer (1994) by using dark hues at each corners of the legend to represent categories that are extremes for both variables. A light and white color was placed at the center of the legend to represent the class that contains the critical value or midpoint of both variables. The remaining colors are lighter than the corners and lightness was adjusted in response to critical values within the data ranges of both variables. Since bivariate mapping is easy to increase cognitive overload, there are no more than three classes applied for bivariate maps in this system. In addition, the data classification was applied with appropriate ranges, such as using week to represent the early and late of the growing season. The design of the map animation is complicated and requires several tests before publishing. The important factor is the dynamic visual variables. This research applied three dynamic visual variables for map animation. Duration was defined for each scene of map animation; the rate of change was adjusted with appropriate time interval of each frame to make smooth change, while it has long enough duration for user to capture the change. All maps were arranged in chronological order. The color scheme and the number of classes are also important for designing map animation. The maps in map animation applied the same technique with thematic map but some maps need adjustments and modifications to display appropriate results on animation. In addition, to facilitate an understanding of the process of change, this system includes interactive temporal control and temporal legend. The temporal legend is designed for users to focus on the detail as well as to depict change over time and space. Temporal control is useful for users to review the events they may have missed within the sequence. With these effective techniques, map animation is able to communicate the concept of change for spatial and temporal context of phenological information. 220 7.2.3 Creating Map Service Map service was created to control a map server on how to generate a map when a request is received. Map service is a standard protocol used to request mapping data from a map server. Moreover, it is also determined and designed by considering the needs of users interactively between web browser and map server. The tool used for creating a map service is UMN MapServer. MapServer is an open source development environment for building spatiallyenabled internet applications. It can be run as a CGI program or via MapScript, which supports several programing languages. In this research, MapServer that is used to develop a map service is the MapServer for windows (MS4W) and it run as CGI Program. The fundamental for creating map service is to create a map service configuration file. It is called as a Mapfile and is used for designing the map elements and configuring the behavior of the map and its layer. The Mapfile defines the relationship among objects as well as points MapServer to where data are located and defines how things are to be drawn (The University of Minnesota, 2011). Generally, each instance of map service needs to have its own Mapfile in which some parameters and some metadata entries are mandatory. For IPA, functions of this system were developed by using Python to automatically create the Mapfile. 7.2.4 Creating Map Animation Map animation in this system was created by using the Flash program. Thematic maps with their legend were created in ArcMap program and then were imported to Flash program. It is necessary to design the layout of animation in advance to indicate the location of map components, such as a title, legend, explanation, etc., as well as the interactive functions. In the Flash program, duration, rate of animation change, and order of maps are set for transition of animation. This system designs and adds a temporal legend with indicators to show the time 221 change (e.g., year). Temporal legend is important for map animation to indicate what time period is currently shown, as well as the total duration of times that will be shown in the animation. Additionally, the script in Flash provides a variety options to add interactivity to animation. This system uses action scripting to develop animation control which allows users to compare data of different years and to view the map at users’ preferred pace. The animations were exported from Flash program in swf format and combined to IPA by using HTML. 7.2.5 Establishing and Customizing the Prototype of IPA After creating the map service, an interface and functions in the IPA prototype are built and customized by HTML and JavaScript. The IPA is developed based on Google Maps API to build an interactive map. Using Google Maps API version 2 as a platform is an easy way for users to display, control, and compare between any layers and other information. Generally, Google Maps API does not support MapServer directly; however, it does support various ways of overlaying images on top of the map. In this research, IPA functions were developed by using JavaScript and Python as a connector to communicate between Google Maps API and MapServer. When a user clicks on a map or a tool, this function automatically generates and sends a request to the map server. When a response is returned from the map server, this function also processes the response by retrieving and displaying a map as an image on the Google Map platform. With this platform, only one image map service can be displayed at a time. In addition to the web map service function, other basic and advanced functions are developed and specified in the system in order to create the map elements, control a map, manipulate the data and display specific information. For example, dynamic time-series graph was developed by using Google Chart API. The web browser must be Firefox to handle the communications for requests, responses, and the specific functions. 222 7.3 Capabilities of IPA To understand the capabilities of IPA, accessibility, data, navigation and interactive functionalities are explained as follows: 7.3.1 Accessibility IPA was developed by UMN MapServer, HTML, Javascript, Python, and Flash. The testing results reveal that IPA could improve public’s accessibility and interaction to phenological information. IPA delivers information to users with fewer requirements for the client-side computing environment. However, this system is the most compatible when running on Firefox browser. Other browsers, such as Microsoft Internet Explorer and Google Chrome, can be used to access IPA but some functions might not be supported by these browsers. 7.3.2 Data The data in this system are listed in Table 3.6. There are maps of EVI, phenological parameters, rainfall seasonality, trend, bivariate, data quality assessment, and GHG emissions (on a local scale) in the thematic map layers. The boundary of Indochina region and Lopburi province in Thailand (case study area) are provided in general map layers. The formats of the data are in raster format for image maps and vector format for supporting layers. Additionally, EVI and seven phenological parameters are selected to be shown in map animation. 7.3.3 Functionality This system builds two groups of interactive functions for users to explore maps. The two groups of functions are basic and advanced functions. 223 1) Basic Functions The basic functions consist of map navigation (pan and zoom in/out), map coordinate, and map-based overlay. These functions are developed by using the Google Maps API. Moreover, the overview map is added for users to see a larger spatial extent so that users can get a better understanding of the spatial relationship of the current map and the entire region. Scale bar is included for distance measurement. In addition, this system has the function to adjust transparency of overlay thematic map. 2) Advance Functions a) Dynamic generating layer content: This system provides the function to automatically and dynamically generate the map layer content corresponding to the maps in data source. When files in database change, map layer content in the interface also correspond to those changes and users can select any updated layers. b) Dynamic legend: This function was developed in this system to automatically change the map according to the map layer selection (with the activated layer). For example, every time user selects a map layer, the system not only generates and provides a thematic map, but also automatically creates a map legend corresponding to map content. In addition, this legend is created as an HTML object that can support the dynamic change and define interactive functions. In some complex maps, such as the bivariate maps, the dynamic legend can also provide more information for users to easily interpret the map. c) Dynamic time-series graph: This function was developed based on Google Chart API and JavaScript. It is performed to access and retrieve the pixel-based information from the temporal images in order to draw a time-series graph. When users click on a specific location on the thematic map, a time-series graph is automatically generated and shown in the popup 224 window at that location, and the graph will be changed when users click on other locations. For example, when users click on specific location on the start date of the growing season map, the system will generate a time-series graph of that location to show start date for ten years. The graph with the x-axis representing years and the y-axis representing different values depend on the map selection. Moreover, the interactive function on the dynamic graph was designed to change the graph to other parameters of that location by simply clicking on the graph without moving back to layer selection. This interactive dynamic graph function is very useful for users to explore and compare the patterns of that parameter during ten years. d) Map animation: Map animation is another effective technique used in IPA. EVI and seven phenological parameters were selected to show an animation in Flash environment. Map animation was generated separately from map server but embed into the system. Map animation is performed smoothly in a specific panel but completely overlay on the panel of map extent. Map animation can be reached by clicking on the animation button in the top toolbar and theme of maps is changed by selecting a specific layer on the drop-down list. Interactive animation control and temporal legend are provided for users to depict and explore change over time and space. 225 7.4 Discussion and Conclusions 7.4.1 Achievement: The Concept of Visual Exploration in IPA The main goal of IPA is to provide visual exploration with an easy access and interactive functions for phenological data and related variables in SEA. From the literature reviews (DiBiase et al., 1992; Blok, 2006), it was determined that the visual exploration needs the tools to explore maps, for example, navigation tool. IPA is mainly developed by integrating MapServer and Google Maps API, which consist of several interactive functions of both basic and advanced functions. A GUI with these tools provides researchers a convenient way to conduct phenological and environmental research. In this system, cartographic techniques and basic and advanced functions are provided for users to explore the map. Basic functions include map navigation, map coordinate, map-based overlay, and layer transparency. Advanced functions consist of dynamic generating layer content, dynamic legend, dynamic time-series graph, and map animation. Additionally, cartographic techniques and designs were applied in this system: thematic map, bivariate map, and map animation. Thematic maps display the patterns of change over time and space. Bivariate map can communicate the relation between two variables. These bivariate maps enhance the ability to detect and comprehend important phenomena and major conclusions. In addition, map animation allows users not only to detect changes but also to understand meanings encoded in transitions within the dynamic displays and with interactive functions. Therefore, IPA provides highly interactive graphic displays that facilitate the environmental analysis of the change across space and through time. This process of exploration with interactive and dynamic interfaces enables 226 users to derive meaningful information as well as to construct knowledge in the research process that is the concept of visual exploration. 7.4.2 Insights: Challenges and Limitations IPA was designed and developed for interactive representation of spatio-temporal environmental maps. The GUI and dynamic and interactive functions of IPA provide easy access to the maps with simple controllers for map navigation and exploration. However, there are challenges and limitation in design and development of this system. In case of technical issues in IPA, the system in IPA was mostly created by Open Source software. Although Open Source software provides interoperable application, some components are in individual software or functions that are developed by different groups leading to different versions. One component featuring additional desired functionality may result in a cascade of updating other components, particularly for a newly release version. In addition, these components must be configured, compiled and installed individually, which makes the maintenance of such systems a time consuming process. In terms of cartographic design, the nature of raster images causes the difficulty in visual exploration. The complexity of raster maps with a larger number of pixels could be misleading for interpretation and could increase cognitive overload for users to be confused with the maps. Consequently, this visualization system attempts to find appropriate color scheme for each variable—for example, different color hue for representing the month of the start and the end of the growing season, or lightness for the quantitative difference of the length of the growing season. These techniques might help user easily understand and detect the change. However, some image maps are difficult to design, such as the trend map that have the variation of values. Although the color schemes of these maps were carefully selected to have only five 227 classes, they provide the poor visualization when displayed on a regional scale. However, these maps show a clear display by using navigation tools to explore site by site. This research also found difficulties in designing bivariate map and map animation. Bivariate maps received considerable criticism for their presumed failure to communicate either information about individual distributions or the correlation between them (Olson, 1981). In addition to the design of data classification, color scheme and appropriate variable selection, this research found that the legend is important to explain the meaning of bivariate map. Therefore, this research added more information on the legend of bivariate maps by adding one word at the corner of the legend to explain extreme conditions (Figure 4.16). For example, the word “shift” was added to one corner of the legend indicating that the growing season starts and ends late. In addition, this system provides the additional explanation of four extreme conditions below the legend. These techniques allow users to understand and easily interpret the map. Additionally, to enhance the effectiveness of bivariate map in order to show spatio-temporal changes, mean difference or mean deviation of two maps can be applied. This research performed the mean difference between phenology of each year with the 10-year mean values based on pixel by pixel to indicate both the relationship between two parameters and the relation between different years. With these techniques, users can analyze the data, detect and comprehend important phenomena, and find a major conclusion from these maps. This research attempts to apply cartographic design to facilitate visualization in map animation such as using appropriate number of classes and applying dynamic visual variables and interactive map animations. However, the evaluation has not been addressed for this visualization prototype. The problem may exist without notice, particularly for map animation. For the problem of the map animation, numerous studies have pointed out that 228 information in an animation can be missed especially, when the data in animation are complex (Slocum et al., 2009). This is because complex data could increase cognitive overload. In addition, the problems of “change blindness” and “change blindness blindness” are critical in map animations (Fish, 2011). These two situations might cause users to miss the changes; moreover, they believe that they interpreted the map correctly (Fish et al., 2011). Therefore, users do not review the display, resulting in the underestimation of changes. With these possible problems, the evaluation is important to improve the performance of this prototype. 229 CHAPTER 8 CONCLUSIONS AND FUTURE ENVISIONS 8.1 Major Findings 8.1.1 Spatio-temporal Variation of Vegetation Phenology and Rainfall Seasonality Phenological patterns derived from MODIS products (16-day EVI with 250 m resolution from 2001-2010) and rainy season observations from TRMM data (daily rainfall with 0.25x0.25 spatial resolution from 2001-2010) are proven to be useful tools for monitoring the response of vegetation to rainfall and identifying the sensitivity of various vegetation formations to climate variability. These data can identify the changes, the relationships between phenology and climate variability, as well as the drivers of phenological changes related to climate variations. However, the study of these patterns has not been addressed in SEA, particularly on a regional scale. In addition, the method used in this research minimized the need for ground based measurement and real-time climatic and environmental data which are unavailable in some locations. To explore vegetation conditions in SEA, this research selected important phenological parameters that can emphasize phenological patterns and changes, including the timing of the growing season, magnitude of vegetation, and vegetation production. This research investigated the patterns of phenological changes in both first and second growing seasons. In addition to phenological parameters, EVI value was also selected as an indicator of the overall greenness of vegetation to study vegetation conditions in this region. This study not only explored the spatio-temporal patterns of phenology on a regional scale, but also the local patterns. This is because local analysis provides valuable information and a meaningful conclusion, particularly for drivers of phenological changes. The hotspots and important 230 representative sites were selected to explain the EVI and phenological patterns as well as the driver of phenological changes due to rainfall variability and human management on a local scale. MODIS-derived estimates of EVI and phenology in this research indicated that EVI and phenological patterns varied spatially according to climate variations and human management. The results showed geographically and ecologically homogeneous patterns in the same ecoregion or biome type. Phenology of naturally vegetated areas and agricultural areas showed differences in timing and magnitudes. The phenological dates occurred later in natural vegetation but started earlier in agricultural areas. When considering the length of the growing season, forest areas have longer growing seasons than agricultural areas. The trend of overall regional mean EVI value in SEA from 2001 to 2010 was decreasing, and phenological trends appeared to shift towards a later and slightly longer growing season in the Peninsula. Regional vegetation dynamics over SEA exhibited patterns that were associated with major climate events such as El Niño in 2005. Various parameters showed distinctly different patterns in years of these extreme climate events. Additionally, the change in timing of growing season with variation of amplitude and integral of forest areas suggests the land use change, particularly in the northern and eastern regions of Vietnam and the southern regions of Laos. Phenology in agricultural areas mainly varied according to human management. The decrease in amplitude and integral of forest areas indicated the loss of high-quality forests, while the increase in amplitude and integral of agricultural areas suggested the improvement of agricultural practices. The rainy season in SEA tends to start early and end late. The length of rainy season was slightly longer; however, the amount of rainfall has decreased from 2001 to 2010. 231 Rainfall patterns were changing regionally, with increases in some locations and decreases in others. Annual rainfall has decreased across most of the western and northern regions of the Peninsula but has increased across the eastern region of Thailand. The relationship between phenological changes and rainfall variability is a key to understand the influence of climate change on biophysical characteristics. The regional data described in this research provide essential information for understanding and modeling the effects of climate variability on ecosystem dynamics. The results reveal a significant positive correlation between SGS and SRS overall, although negative correlations can be found in some areas that are used most extensively for agriculture. The forest areas and rainfed croplands show similar patterns, in which the growing season shifted according to changes in rainfall, especially in Myanmar. The rate of phenological changes with respect to rainfall seasonality varied among different land cover types and ecosystems. Although the trends of phenology were different from place to place, the major drivers of the phenological changes are climate and human management. The regional and local scale analyses indicate that rainfall is a dominant force in naturally vegetated areas and rainfed croplands, whereas human management is a key factor in heavily agricultural areas with irrigated systems. However, it was apparent that some croplands (e.g., sugar cane, cassava) also need rainfall to control vegetation growth. The phenology and rainfall were highly correlated in these places. In order to quantify the agreement between the satellites derived phenology and field observations as well as satellite TRMM rainfall and station rainfall, this research used the data from field survey in Thailand to compare with MODIS phenology and applied rainfall station in Thailand to compare with TRMM rainfall. The comparison demonstrated the 232 acceptable agreement between satellite data and field data; although discrepancies exist in some cases of satellite phenology and field observations due to the different technique used to acquire phenology. In summary, the general pattern and phenological changes in SEA at a regional scale apparently exhibited relative correspondence with rainfall seasonality and major climate events, particularly in El Niño and La Niña years. Climate variability has affected SEA agriculture, leading to greater losses in economics of this region as well as posing a serious threat to future food security. Not only has climate variability affected phenological changes, human management, with land use practices, can also profoundly affect land surface phenology. Changes in land use practices may contribute to large-scale alterations in the timing of landsurface/atmosphere boundary conditions through phenological changes. Additionally, the results of this research has identified the hot spots or sensitive areas of decreasing EVI and forest loss as well as significant changes of phenology in croplands. These results could guide management priorities for protection and conservation in forest areas and agricultural management in cropland areas. As SEA is more vulnerable to the impacts of climate change, particularly in agricultural areas, this meaningful information is useful for monitoring and assessing the effects of climate change on ecosystem dynamics and environmental changes at both regional and local scales. 8.1.2 Spatio-temporal patterns of GHG Emissions from the DNDC Model under Different Scenarios This research demonstrates the potential integration of remote sensing mapping with the biogeochemical model to quantify spatio-temporal patterns of GHG emissions (CH4, CO2, and N2O) from rice fields in Thailand for both single and double rice cropping systems. This research took a new step to advance the regional application of DNDC by using regional 233 database based on a grid-based system instead of average input parameters in a large spatial unit (e.g., county unit). Thus, the simulation of grid based system at the site mode of DNDC can differentiate the difference of soil properties, climate data, and farming management and can substantially improve the accuracy of the GHG estimations from the DNDC model. Additionally, phenological information is indispensable for DNDC but it is difficult to access or is often unavailable. Satellite time-series data can provide such important information. This research used spatial and temporal characteristics of phenology derived from remote sensing data for DNDC to quantify emissions. The results of model simulation indicate that the new database system and spatio-temporal data can provide more accurate emission predictions from the DNDC model. Spatial and temporal dynamics of GHG emissions were characterized under different scenarios to investigate the effects of the various factors on GHG emissions from rice fields, including climate and human management (fertilization and phenology). The results indicate that paddy rice fields were a source of atmospheric GHG. The change of GHG emissions was influenced by human management, climate variation, and physical geography. The findings reveal that phenology is the main factor affecting the changes in CH4 emissions. The change of CO2 emissions was relatively smaller than CH4 and N2O in all scenarios but showed rapid changes in extreme climate years. N2O emissions were strongly related to climate variation, especially rainfall changes. Soil texture plays an important role on CH4 emissions with a high level of emissions in light texture soil. Rice intensification with longer length of flooding period and high application rates of fertilizer extremely induce CH4 production due to the change in the chemical 234 and biophysical properties of the plant-soil system. Additionally, standard emission factors of CH4 emissions obtained from the IPCC approach and research in Thailand were used to compare with DNDC results. The modeled emission results were in good agreement with these standard emission factors. In conclusion, human management in farming practices apparently affects spatiotemporal patterns of GHG emissions. The results of this research suggest that practical mitigation options should be carefully regulated to balance the emission of different gas types more efficiently, as well as to maintain and improve grain yields. This is because a decrease in one gas emission may lead to significant increases in other gases’ emissions. The results imply that mitigation of CH4 emissions should be emphasized for rice fields because its emission is much higher than CO2 and N2O. These gas emissions not only represent a negative impact on the environmental quality but could also lead to economic losses. Balancing between optimum crop yield and environmental safety is a great challenge for agricultural management. Integration of remote sensing and a biogeochemical model provides the opportunity to identify mitigation strategies to balance food production and environmental protection as well as to predict the impacts of climate change on agroecosystem. 8.1.3 The Interactive Phenological Atlas for SEA This research developed improved visualization methods for phenological information dissemination by using the concepts of visual exploration and Internet atlas. The prototype of the Interactive Phenological Atlas for SEA (IPA) was established to display thematic map, bivariate map, and map animation using web mapping application. Phenological information and related environmental variables (EVI, phenological parameters, rainfall 235 seasonality, and annual data quality assessment maps in SEA) were processed to represent spatial and temporal patterns of environmental changes in this region. Additionally, spatial distributions and changes of GHG emissions on a local scale analysis were also included in this system. The objective of this system is data exploration. This system allows users to derive meaningful information from interactive processes. The target users are researchers in environmental and remote sensing fields. The functions in this system enable users to navigate the maps in both spatial and temporal extents. The basic functions provide tools for map navigation, map coordinate, map-based overlay, and layer transparency, while the advanced functions include dynamic generating layer content, dynamic legend, dynamic time-series graph, and map animation. Cartographic design was carefully applied to each map according to its type (thematic maps, bivariate map, and map animation). This research selected appropriate techniques and designs in order to effectively communicate complex spatio-temporal geographic phenomena, such as visual variables, data classification, and number of classes for thematic mapping, and dynamic visual variables for map animation. The most important outcome of this system is that it allows users to explore the complexity of dataset. This system facilitates users with an interactive graphic user interface and effective functions in order to provide a convenient way for researchers to conduct phenological research. The visualization system for this application has little availability and little exploration in the field. Consequently, this phenology and climate variation exploration could generate research ideas and confirm information in research process that would be useful for environmental study. 236 8.2 Challenges and Limitations Validation is a significant concern in remote sensing-based analysis. The method of validation is very challenging due to the scale differences, time consumption, and high cost. Although comparisons of MODIS-derived estimates of phenology and ground observations in this research provide satisfactory results, substantially more field data are required to fully assess the phenology based remote-sensing time-series data. Additionally, field observations of phenology in a variety of environments are required and field data collection should be specifically designed to validate remote sensing phenology. Satellite derived phenology is useful for the environmental study. However, their applications are hindered by atmospheric conditions, particularly in the tropical zone, due to frequent cloud cover in the wet season (growing season) and fire-induced aerosols in the dry season. A filtering method is needed for time-series data before extracting phenology. However, this approach is sometimes insufficient to eliminate cloud contamination resulting in unsuccessful phenological extraction in some pixels. Therefore, annual data quality assessment maps are provided in this research for comparing the level of quality associated with annual EVI and annual phenological patterns. The new MODIS EVI collection 6 will be released in the near future, which could provide more accurate data that would be useful for phenology research. Although MODIS satellite time-series data with high temporal resolution is free, the 250 m resolution leads to the spatial heterogeneity of land cover classes in one pixel. These mixed pixels introduce error into the analysis of interannual variability. The interpretation of remotely sensed phenological parameters at this spatial resolution should be addressed. The challenge of GHG emissions quantifying from DNDC model is the database development. Due to numerous input data for DNDC model, it is required to use default values 237 provided in the DNDC model or estimated values to run the model because some data are unavailable or lack of information from field observations or data records. This may lead to incorrect results in some cases. Additionally, the validation is also challenging for the DNDC modeling. The validation with field data is necessary for the model performance but there is no field data available in Thailand. IPCC approach and standard emission factors were applied in this research to compare with modeled results. However, the standard emission factor is unable to reflect the change of input parameters (e.g., climate, soil, management practices) resulting in the discrepancy between the simulated results and the standard factor. Both technical issues and cartographic designs are challenging for IPA. The Open Source software applied in this system provides a great opportunity to develop and modify the functions but it is difficult to keep up-to-date because of their various different versions. Additionally, cartographic design should be carefully considered for effective visualization. Different techniques and several tests are required to design thematic map, bivariate map, and map animation for visualizations. To improve the performance of this visualization prototype, an evaluation is important to demonstrate the success of the system. Evaluation can indicate the advantages and drawbacks of functions as well as the performance of the system in order to develop effective systems for users. 8.3 Future Research The data utilized in this study were analyzed for a ten-year period, which may not be long enough to assess long-term ecosystem changes. Future ecological and environmental change analyses could benefit from extended study periods or integration with previous research to assess ecosystem dynamics as well as global climate change. 238 Future work is needed to validate remote sensing-based observations with high-resolution imagery or field observations. A set of field observations collected with effective techniques should be developed to provide a better basis for the comparison between field measurements and remotely sensed observations. This assessment is useful for the satellite-derived phenology research. Influence of soil properties and other climate data should be considered in future research. Although rainfall is a significant factor for phenological changes in this region, vegetation green-up is also associated with soil properties, particularly soil moisture, which is related to rainfall and vegetation types. Temperature could alter changes in growing season in some locations. The study of these factors may contribute to the understanding of the ecosystem dynamics and their disturbances for future adaptation strategy development in this region. Disaster is the critical issues in this region and it requires the monitoring system. Phenology based on remote sensing data could be used for estimating, indicating, and monitoring vegetation condition for environmental disasters (e.g., drought, flood, forest fire), especially in sensitive areas. This analysis could contribute to ecological effects of environmental change and the impact of climate and human on the environment due to stressful conditions. Integration of DNDC and remote sensing technique used in this study can be applied to other rice areas in this region over large-scale areas at country or regional levels because of lower cost and less time consuming. Additionally, future research should conduct sensitivity analysis using different inputs or major drivers, such as soil, to explore variations in GHG emissions from rice fields. These studies would be useful and highly valuable for regional and global emission inventory. 239 Future improvements of the IPA prototype should primarily address the performance for rendering image maps. Additional effective techniques for map animation could enable visual exploration. For example, visual salience can be applied on the complex maps by highlighting or signaling important content with flashing or dominant colors. Additionally, the synchronization of map animations allows identification of pattern correspondences in time-series by displaying more than one animation on the screen. This technique would contribute to the understanding of the relationship between related parameters in the spatial and temporal contexts. 240 APPENDICES 241 APPENDIX A EVI and Phenological Parameters 2001-2010 Figure A1 EVI and phenological parameters 2001-2010 242 Figure A2 Start date of the Growing Season 2001-2010 (first growing season) 243 Figure A3 End date of the Growing Season 2001-2010 (first growing season) 244 Figure A4 Length of the Growing Season 2001-2010 (first growing season) 245 APPENDIX B Rainfall Seasonality 2001-2010 Figure B1 Mean Annual Rainfall 2001-2010 246 Figure B2 Start Date of the Rainy Season 2001-2010 247 Figure B3 End Date of the Rainy Season 2001-2010 248 Figure B4 Length of the Rainy Season 2001-2010 249 APPENDIX C Data Quality Assessment Map Figure C1 MODIS EVI data quality Assessment 2001-2010 250 Figure C2 MODIS EVI data quality assessment (10-year mean) 251 APPENDIX D Multi-year Simulation of GHG Emissions from DNDC Model (Scenario 3-5) Figure D1 CH4 emission rate in scenario 3 252 Figure D2 CO2 emission rate in scenario 3 253 Figure D3 N2O emission rate in scenario 3 254 Figure D4 CH4 emission rate in scenario 4 255 Figure D5 CO2 emission rate in scenario 4 256 Figure D6 N2O emission rate in scenario 4 257 Figure D7 CH4 emission rate in scenario 5 258 Figure D8 CO2 emission rate in scenario 5 259 Figure D9 N2O emission rate in scenario 5 260 REFERENCES 261 REFERENCES Aden, C., Schmidt, G., Schonrock, S, & Schroder, W. (2010). Data analyses with the WebGIS WaldIS. European Journal of Forest Research, 129, 489-497. Aggarwal, P. K., Kalra, N., Chander, S., & Pathak, H. (2004). InfoCrop: A generic simulation model for annual crops in tropical environments. Indian Agricultural Research Institute, New Delhi. ArcGIS Resource Center. (2011). Thiessen polygons. Retrieved from http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00080000001m000000) Asian Development Bank (ADB). (2009). Climate Change in Southeast Asia: A Regional Review. Manila, Philippines: Asian Development Bank. Bachoo, A & Archibald, S. (2007). Influence of using date-specific values when extracting phenological metrics from 8-day composite NDVI data. International workshop on the analysis of multi-temporal remote sensing images, 2007. Multitemp 2007 IEEE Conference. Leuven, Belgium, 18-20 July 2007, pp 4. Beck, P.S.A., Atzberger, C., & Hogda, K.A. (2006). Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI. Remote Sensing of Environment, 100(3), 321-334. Beck, P.S.A., Jönsson, P., Hogda, K.A., Karlsen, S.R., Eklundh, L., & Skidmore, A.K. (2007). A ground-validated NDVI dataset for monitoring vegetation dynamics and mapping phenology in Fennoscandia and the Kola peninsula. International Journal of Remote Sensing, 28(19), 4311-4330. Blok, C. A. (2006). “Interactive animation to visually explore time series of satellite imagery.” Book Series: Lecture Notes in Computer Science: V3736: 71-82, Springer-Verlag Berlin. Borchert, A. (1999). Multimedia atlas concepts. In Cartwright, W., Peterson, M.P., &Gartner, G. (Eds), Multimedia cartography (pp. 75-86). Germany: Springer. Boonpragob, K. 2005. Crisis or Opportunity: Climate Change Impacts and Thailand. Greenpeace Southeast Asia, Thailand. Bradley, B.A., Jacob, R.W., Hermance, J.F., & Mustard, J.F. (2007). A curve fitting procedure to derive inter-annual phonologies from time series of noisy satellite NDVI data. Remote Sensing of Environment, 106, 137-145. 262 Bradley B.A. & Mustard J.F. (2008). Comparison of phenology trends by land cover class: a case study in the Great Basin, USA. Global Change Biology, 14, 334-346. Brewer, C.A. (1994). Color use guidelines for mapping and visualization. In A.M. MacEachren and D.R.F. Taylor (Ed.), Visualization in modern cartography (pp. 123-147). New York: Pergamon. Bridhikitti, A., & Overcamp, T.J. (2012). Estimation of Southeast Asian rice paddy areas with different ecosystems from moderate-resolution satellite imagery. Agricuture, Ecosystems & Environment, 146(1), 113-120. Brown, L., Javis, S. C., & Headon, D. (2001). A farm-scale basis for predicting nitrous oxide emissions from dairy farms, Nutrient Cycling in Agroecosystems, 60, 149– 158. Butterbach-Bahl, K., Kesik, M., Miehle, P., Papen, H., & Li, C. (2004). Quantifying the regional source strength of N-trace gases across agricultural and forest ecosystems with process based models. Plant and Soil, 260(1/2), 311-329. Cai, Z, Sawamoto,T, Li, C., Kang, G., boonjawat, J., Mosier, A., Wassmann, R., & Tsuruta, H. (2003). Field validation of the DNDC model for greenhouse gas emissions in East Asian cropping systems. Global Biogeochemical Cycles, 17 (4), 18-1-18-10. Carrao, H., Goncalves, P., & Caetano, M. (2010). A nonlinear harmonic model for fitting satellite image time series: analysis and prediction of land cover dynamics. IEEE Transactions on Geoscience and Remote Sensing, 48(4), 1919-1930. Cao, M., Dent, J. B., & O. W. Heal. (1995). Modeling of methane emission from rice paddies. Global Biogeochemical Cycles, 9, 183– 195. Chen, J., Jonsson, P., Tamura, M., Gu, Z., Matsushita, B., & Eklundh, L. (2004). A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sensing of Environment, 91, 332-344. Choi, J.Y., Engel, B.A., Theller, L., & Harbor, J. (2005). Utilizing web-based GIS and SDSS for hydrological land use change impact assessment. American Society of Agricultural and Biological Engineers, 48(2), 818-822. Chokngamwong, Roongroj & Chiu, L.S. (2008). Thailand daily rainfall and comparison with TRMM products. Journal of Hydrometeorology, 9(2), 256-266. 263 Cleland, E.E., Chuine, I., Menzel, A., Mooney, H.A., & Schwartz, M.D. (2007). Shifting plant phenology in response to global change. Trends in Ecology and Evolution, 22(7), 357365. Cooke, R. & Toda, A. (2008). Confronting Climate Change and Land Degradation in Viet Nam: Increasing Finance for Sustainable Land Management. A collaborative endeavour of the Global Mechanism of the United Nations Convention to Combat Desertification and the International Fund for Agricultural Development. Retrieved from http://www.globalmechanism.org/dynamic/documents/... Corlett, R.T. & Lafrankie, J.V. (1998). Potential impacts of climate change on tropical Asian forests through an influence on phenology. Climatic Change, 39, 439-453. Cox, D.J. (1990). The art of scientific visualization. Academic Computing 4(6), 20-22. Cui, Y.P., Liu, J.Y., Hu, Y.F., Kuang, W.H., & Xie, Z.L. (2012). An analysis of temporal evolution of ndvi in various vegetation-climate regions in Inner Mongolia, China. Environmental Sciences, 13, 1989 – 1996. de Beurs, K. M. & Henebry, G. M. (2004). Land surface phenology, climatic variation, and institutional change: Analyzing agricultural land cover change in kazakhstan. Remote Sensing of Environment, 89(4), 497-509. de Beurs, K.M. & Henebry, G.M. (2005). A statistical framework for the analysis of long image time series. International Journal of Remote Sensing, 26(8), 1551-1573. de Beurs, K.M. , Wright, C.K., & Henebry, G.M. (2009). Dual scale trend analysis for evaluating climatic and anthropogenic effects on the vegetated land surface in russia and kazakhstan. Environmental Research Letters, 4(4), 045012. de Beurs, K. M., & Henebry, G. M. (2010). Spatio-temporal statistical methods for modelling land surface phenology. In Hudson, I. L. & Keatley, M. R (Eds.) Phenological research (pp. 177-208): Springer Netherlands. Del Grosso, S., Ojima, D., Parton, W., Mosier, A., Peterson, G., & Schimel, D. (2002). Simulated effects of dryland cropping intensification on soil organic matter and greenhouse gas exchanges using the DAYCENT ecosystem model. Environmental Pollution, 116, S75– S83. DiBiase, D., MacEachren, A., Krygier, J., & Reeves, C. (1992). Animation and the Role of Map Design in Scientific Visualization. Cartography and Geographic Information Systems 19(4), 201-214. 264 Didan, K. & Huete, A. (2006). MODIS Vegetation Index Product Series Collection 5 Change Summary. Retrieved from http://modis-250m.nascom.nasa.gov/QA_WWW/forPage/ MOD13_VI_C5_Changes_Document_06_28_06.pdf Eiumnoh, A., Jugsujinda, A., Parkpian, P., Shrestha, R.P., Baimoung, S., Pongsai, S., Kumnerdmanee, S., Pattanapong, B., Tanchot, T., Tongkailad, R., & Poonsawad, V. (2002). Development of Database for Greenhouse Gases Emission with Reference to Methane in Thailand. The Thailand Research Fund (TRF): research report. EPA (United States Environmental Protection Agency). (2006). Data quality assessment: statistical methods for practitioners. Washington, DC: Author. Retrieved from http://www.epa.gov/QUALITY/qs-docs/g9s-final.pdf Fish, C., Goldsberry, K.P., & Battersby, S.E. (2011). Change blindness in animated choropleth maps: an empirical study. Cartography and Geographic Information Science, 38(4), 350362. Fisher, J. I., & Mustard, J. F. (2007). Cross-scalar satellite phenology from ground, Landsat, and MODIS data. Remote Sensing of Environment, 109(3), 261-273 Food and Agriculture Organization of the United Nations (FAO). (2007). Adaptation to climate change in agriculture, forestry and fisheries: Perspective, framework and priorities. Rome, Italy: Food and Agriculture Organization of the United Nations. Fumoto, T., Kobayashi, K., Li, C., Yagi, K., & Hasegawa, T. (2007). Revising a process-based biogeochemistry model (DNDC) to simulate methane emission from rice paddy fields under various residue management and fertilizer regimes. Global Change Biology, 14(2), 382-402. Fumoto, T., Kobayashi, K., Li, C., Yagi, K., & Hasegawa, T. (2008). Revising a process-based biogeochemistry model (DNDC) to simulate methane emission from rice paddy fields under various residue management and fertilizer regimes. Global Change Biology, 14, 382–402. Gale, G.A., Towprayoon, S., & Smakgahn, K. (2005). Development of a database for estimating methane emissions from rice fields in Thailand. The Thailand Research Fund (TRF): research report. Gartner, G. (1999). Multimedia GIS and the web. In Cartwright, W., Peterson, M.P., &Gartner, G. (Eds), Multimedia cartography (pp. 305-314). Germany: Springer. 265 Goldsberry, K.P. & Battersby, S. (2009). Issues of change detection in animated choropleth maps. Cartographica 44(3), 201-215. Hachler, T. (2003). Online visualization of spatial data a prototype of an open source Internet map server with back end spatial database for The SWISS National Park (Master Thesis). Retrieved from http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.138.8561 Harrower, M. (2002). Visualizing Change: Using Cartographic Animation to Explore RemotelySensed Data. Cartographic Perspective 39, 30-42. Harrower, M. (2008). The Role of Map Animation for Geographic Visualization. In Martin Dodge, Mary McDerby, & Martin Turner (Eds.), Geographic visualization (pp. 49-65). NJ: John Wiley & Sons, Ltd. Hermance, J. F., Jacob, R. W., Bradley, B. A., & Mustard, J. F. (2007). Extracting phenological signals from multiyear AVHRR NDVI time series: Framework for applying high-order annual splines with roughness damping. IEEE Transactions on Geoscience and Remote Sensing, 45(10), 3264-3276. Hirsch, R. M., Slack, J. R.,& R. A. Smith (1982), Techniques of trend analysis for monthly water quality data, Water Resources Research, 18(1), 107–121. Hirsch, R. M. & Slack, J. R. (1984). A nonparametric trend test for seasonal data with serial dependence. Water Resources Research, 20(6), 727-732. Hird, J. N., & McDermid, G. J. (2009). Noise reduction of ndvi time series: An empirical comparison of selected techniques. Remote Sensing of Environment, 113(1), 248. Horwath, W.R. (2011). Greenhouse gas emissions from rice cropping systems. In Guo, L., Gunasekara, A.,S., & McConnell, L., L (Eds.), Understanding Greenhouse Gas Emission from Agricultural Management (pp. 67-89). American Chemical Society: Washington, DC. Huete, A., Didan, K., Miura, T., Rodriguez, E.P., Gao, X., & Ferreira, L. G. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83(1-2), 195-213. Huffman, G.J., Adler, R.F., Bolvin, D.T., Gu, G., Nelkin, E.J., Bowman, K.P., Hong, Y., Stocke, E.F., Wolf, D.B. (2007). The TRMM Multi-satellite Precipitation Analysis: Quasiglobal,multi-year, combined-sensor precipitation estimates at fine scale. Journal of Hydrometeorology, 8(1), 38-55. 266 IPCC. (1996). Revised 1996 IPCC Guildelines for National Greenhouse Gas Inventories Workbook, Vol. 2, Inter governmental Panel on Climate Change, Cambridge University Press, Cambridge, 1997. IPCC. (2006), 2006 IPCC Guidelines for National Greenhouse Gas Inventories (Volume 4 Agriculture, forestry and Other Land Use), Prepared by the National Greenhouse Gas Inventories Programme, Eggleston H.S., Buendia L., Miwa K., Ngara T. and Tanabe K. (eds). IGES, Japan. IPCC. (2007). Climate Change 2007, Contribution of Working Groups I and II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. United Kingdom and New York, NY, USA: Cambridge University Press Cambridge. Jagadeesh Babu, Y., Li, C., Frolking, S., Nayak, D. R., Datta, A., & Adhya, T. K. (2005). Modelling of methane emissions from rice-based production systems in india with the denitrification and decomposition model: Field validation and sensitivity analysis. Current Science, 89(11), 1904-1912. Jagadeesh Babu, Y., Li, C., Frolking, S., Nayak, D.R., & Adhya, T.K.. (2006). Field validation of DNDC model for methane and nitrous oxide emissions from rice-based production systems of India. Nutrient Cycling in Agroecosystems, 74, 157-174. Jermsawatdipong, P., Vongmaneeroj, A., Kannuj, L., Chairoj, P., Khunathai, H., & Charoensilp, N. (2002). Methane emission from paddy fields in Thailand and data base management. The Thailand Research Fund (TRF): research report. Jong, R.D., Bruin, S., W, A., Schaepman, M.E., & Dent, D.L. (2011).Analysis of monotonic greening and browning trends from global NDVI time-series. Remote Sensing of Environment, 115, 692–702. Jönsson, P., & Eklundh, L. (2004). TIMESAT - a program for analyzing time-series of satellite sensor data. Computers & Geosciences, 30(8), 833-845. Jönsson, P., & Eklundh, L. (2010). TIMESAT 3.0 Software Manual. Department of Earth and Ecosystem Sciences, Lund University and Center for Technology Studies, Malmö University: Sweden. Julien, Y., & Sobrino, J. A. (2009). Global land surface phenology trends from GIMMS database. International Journal of Remote Sensing, 30(13), 3495-3513. Julien, Y., Sobrino, J.A., Mattar, C., Ruescas, A.B., Jimenez-munoz, J.C., Soria, G., et al. (2011). Temporal analysis of normalized difference vegetation index (NDVI) and land surface 267 temperature (LST) parameters to detect changes in the Iberian land cover between 1981 and 2001. International Journal of Remote Sensing, 32(7), 2057-2068. Jutakorn, J. (2011). Onset, Ending, Length and Rainfall Variation for Agriculture in Thailand. Agrometeorological Sector, Meteorological Development Bureau, Meteorological Department, Thailand: research report. Kang, S., Running S. W., Lim, Jong-Hwan., Zhao, M., Park, Chan-Ryul., & Loehman, R. (2003). A regional phenology model for detecting onset of greenness in temperate mixed forests, Korea: an application of MODIS leaf area index. Remote Sensing of Environment, 86(2), 232-242. Keereerom, K. (2011). Estimateing Methane Emission by Using Vegetation Indices from LANDSAT-5 TM : A Case Study of Amphoe Bang Num Priao, (MSc thesis), National Institute of Development Adminsitration, Thailand. Kramer, K., Leinonen, I., & Loustau, D. (2000). The importance of phenology for the evaluation of impact of climate change on growth of boreal, temperate and Mediterranean forest ecosystems: An overview. International Journal of Biometeorology, 44(2), 67-65. Latifovic, R., & D. Pouliot, D. (2007). Analysis of climate change impacts on lake ice phenology in Canada using the historical satellite data record. Remote Sensing of Environment, 106 (4),492-507. Lasco, R.D., Habito, C.M., Delfino, R.J., Pulhin, F., & Concepcion, R.N. (2011). Climate change adaptation for smallholder farmers in Southeast Asia. World Agroforestry Centre, Philippines. Retrieved from http://www.seachangecop.org/files/documents/Lasco_et_al_2011_CCA_Guidebook_v5.p df Lee, R., Yu, F., Price, K.P., Ellis, J., & Shi, P. (2002). Evaluating vegetation phenological patterns in Inner Mongolia using NDVI time-series analysis. International Journal of Remote Sensing, 23(12), 2505-2512. Li, C., Frolking, S., & Frolking T. A. (1992). A model of nitrous oxide evolution from soil driven by rainfall events: 1. Model structure and sensitivity. Journal of Geophysical Research, 97, 9759– 9776. Li, C., Frolking, S., & Harriss, R. (1994). Modeling carbon biogeochemistry in agricultural soils. Global Biogeochemical Cycles, 8(3), 237–254. Li, C., Narayanan, V., & Harriss, R. (1996). Model estimates of nitrous oxide emissions from agricultural lands in the United States. Global Biogeochem. Cycles, 10, 297–306. 268 Li, C. (2000). Modeling trace gas emissions from agricultural ecosystems. Nutrient Cycling in Agroecosystems, 58, 259-276. Li, C., Zhuang,Y. H., Cao, M. Q., Crill, P. M., Dai, Z. H., Frolking, S., Moore, B., Salas, W., Song, W. Z., & Wang, X. K. (2001). Comparing a national inventory of N2O emissions from arable lands in China developed with a process-based agro-ecosystem model to the IPCC methodology. Nutrient Cycling in Agroecosystems, 60, 159–175. Li, C., Zhuang, Y, Frolking, S., Galloway, J., Harriss, R., Moore, B., Schimel, D., & Wang, X. (2003). Modeling soil organic carbon change in croplands of China. Ecological Applications, 13(2), 327-336. Li, C., Mosier, A., Wassmann, R., Cai, Z., Zheng, X., Huang, Y., Tsuruta, H., Boonjawat, J, & Lantin, R. (2004). Modeling greenhouse gas emissions from rice-based production systems: sensitivity and upscaling. Global Biogeochemical Cycles, 18, 1-19. Li, C., Frolking, S., & Butterbach-bahl, K. (2005a). Carbon sequestration in arable soils is likely to increase nitrous oxide emissions, offsetting reductions in climate radiative forcing. Climate Change, 72, 321-338. Li, C., Frolking, S., Xiao, X., Moore, B., Boles, S., Qiu, J., Huang, Y., Salas, W., & Sass, R. (2005b). Modeling impacts of farming management alternatives on CO2, CH4, and N2O emissions: A case study for water management of rice agriculture of China. Global Biogeochemical Cycles, 19, 1-10. Li, Z. & Xiao-Ling, C. (2007). spati-temporal changes ndvi relations precipitation. In Proceedings of Geocomputation 2007, Ireland. LP DAAC (Land Processes Distributed Active Archive Center). (2010). In Vegetation Indices 16-Day L3 Global 250m. Retrieve September 25, 2010, from USGS Online: https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/vegetation_indices/16_day _l3_global_250m/myd13q1 Luus, K.A., 2009. Analyzing pan-Arctic 1982-2006 trends in temperature and bioclimatological indicators (productivity, phenology and vegetation indices) using remote sensing, model and field data. MS thesis. University of Waterloo. Majumdar, D. (2003). Methane and nitrous oxide emission from irrigated rice field: proposed mitigation strategies. Current Science, 84(10), 1317-1326. 269 Mall, R. K. & Aggarwal, P. K. (2002). Climate change and rice yields in diverse agroenvironments of India. I. Evaluation of impact assessment models. Climatic Change, 52, 315–330. Martínez, B., & Gilabert, M. A. (2009). Vegetation dynamics from NDVI time series analysis using the wavelet transform. Remote Sensing of Environment, 113(9), 1823-1842. Mathiyalagan, V., Grunwald, S., Reddy, K.R., & Bloom, S.A. (2005). A WebGIS and geodatabase for Florida’s wetlands. Computers and Electronics in Agriculture, 47, 69-75. Matthews, R. B., Wassmann, R., & Arah, J. (2000). Using a crop/soil simulation model and GIS techniques to assess methane emissions from rice fields in Asia: I. Model development, Nutrient Cycling in Agroecosystems, 58, 141–159. Meng, M., Ni, J., & Zong, M. (2011). Impacts of changes in climate variability on regional vegetation in China: NDVI-based analysis from 1982 to 2000. Ecological Research, 26, 421-428. Minnesota State University & Mankato Water Resources Center and Minnesota Pollution Control Agency. (2009). Minnesota River Basin Statistical Trend Analysis. USA: Retrieved from http://mrbdc.wrc.mnsu.edu/reports/statistical_trends/pdfs/statistical_trends.pdf Moulin, S., Kergoat, L., Viovy, N., & Dedieu, G. (1997). Global-scale assessment of vegetation phenology using NOAA/AVHRR satellite measurements. Journal of Climate, 10(6), 1154-1170. Myneni, R. B., Keeling, C. D., Tucker, C. J., Asrar, G., & Nemani, R. R. (1997). Increased plant growth in the northern high latitudes from 1981 to 1991. Nature, 386(6626), 698-702. National Intelligence Council (NIC). (2009). Southeast Asia and pacific islands: the impact of climate change to 2030: a commissioned research report. Retrieved from http://www.dni.gov/nic/PDF_GIF_otherprod/climate_change/cr201002_southeast_asia_p acific_islands_climate_change.pdf Ninh, N.H., & Kelly, M. (2000). The Impact of El Niño and La Niña On Southeast Asia: the Human Dimension, Policy Lessons and Implications for Global Change. Final Report An Indochina Global Change Network project funded by the Asia-Pacific Network for Global Change Research. Retrieved from http://www.tiempocyberclimate.org/annex/igcn/igcn2000.htm 270 NASA. (2011). TRMM: Tropical Rainfall Measuring Mission. Retrieved from http://trmm.gsfc.nasa.gov/ Ngwenya, F. (2006). Water quality trends in the Eerste River,Western Cape, 1990–2005. MSc thesis, University of the Western Cape, South Africa. Nicholson, S.E., Davenport, M.L., & Malo, ADA.R. (1990). A comparison of the vegetation response to rainfall in the Sahel and East Africa, using normalized difference vegetation index from NOAA AVHRR. Climatic Change, 17, 209-241. Office of Agricultural Economices. (2010). Agricultural Statistics of Thailand. Annual Report. Ministry of Agriculture and Cooperatives, Thailand. Olander, L.P., Malin, D. (2010). Comparison of three biogeochemical process models for quantifying greenhouse gas effects of agricultural management. USA. Retrieved from http://sustainablefood.org/images/stories/pdf/Comparison_of_Three_Biogeochemical_ Process_Models_report.2a.pdf Olson, J.M. (1981). Spectrally encoded two-variable maps. Annuals Association of American Geographers, 71(2), 259-276. Parton, W.J., Stewart, J.W.B., & Cole, C.V. (1988). Dynamics of C,N,P and S in grass soils: a model. Biogechemistry, 5, 109-131. Pathak, H., Li, C., & Wassmann, R. (2005). Greenhouse gas emissions from Indian rice fields: calibration and upscaling using the DNDC model. Biogeosciences, 2, 113-123. Qin, C. (2011). Assessing Phenological Changes and Drivers in East Africa From 1982 To 2006, (MSc thesis), Michigan State University. Reed, B. C., Brown, J. F., VanderZee, D., Loveland, T. R., Merchant, J. W., & Ohlen, D. O. (1994). Measuring phenological variability from satellite imagery. Journal of Vegetation Science, 5(5), 703-714. Reed, B. C., White, M., & Brown J. F. (2003). Remote sensing phenology. In Schwartz, Mark D. (Ed.), Phenology: an integrative environmental science (pp. 363-382). The Netherlands: Kluwer Academic Publishers. Reed, B. C., Schwartz, M. D., & Xiao, X. (2009). Remote sensing phenology: Status and the way forward.. In Noormets, A. (Ed.), Phenology of ecosystem processes (pp. 231-245). New York: Springer. 271 Richard, D. (1999). Web Atlases- Internet Atlas of Switzerland. In Cartwright, W., Peterson, M.P., &Gartner, G. (Eds), Multimedia cartography (pp. 113-118). Germany: Springer. Robeson S.M. (2004). Trends in time-varying percentiles of daily minimum and maximum temperature over North America. Geophysical Research Letters, 31, 1-4. Sakamoto, T., Van Nguyen, N., Ohno, H., Ishitsuka, N., & Yokozawa, M. (2006). Spatiotemporal distribution of rice phenology and cropping systems in the mekong delta with special reference to the seasonal water flow of the mekong and bassac rivers. Remote Sensing Of Environment, 100(1), 1-16. Salas, William, Li, Changsheng, Mitloehner, Frank, and John Pisano. (2008). Developing and Applying Process‐based Models for Estimating Greenhouse Gas and Air Emission From California Dairies. California Energy Commission, PIER Energy‐ RelatedEnvironmental Research. Retrieved from http://www.energy.ca.gov/2008publications/CEC-500-2008-093/CEC-500-2008093.PDF Schwartz, Mark D. (2003). Introduction. In Schwartz, Mark D. (Ed.), Phenology: an integrative environmental science (pp. 3-8). The Netherlands: Kluwer Academic Publishers. Schwartz, M.D. (2006). Onset of spring starting earlier across the Northern Hemisphere. Global Change Biology, 12, 343-351. Siriratririya, O., Patanaponpaiboon, P., Kupkanchanakul, T., Zungsontiporn, S., Iamsupasit, S., Nenyod, B., & Karnjanaphun, S. (2001). Impact of mitigation methane emission from rice field on rice cultivation and productivity in Thailand. The Thailand Research Fund (TRF): research report. Smakgah, K. (1999). Estimating methane emissions from rice fields using emission factors and a Geographical Information System, (MSc thesis), The Joint Graduate School of Energy and Environment at King Mongkut's University of Technology Thonburi, Thailand. Smakgah, K. (2003). Empirical model and process base model application on methane and nitrous oxide emissions from Thai rice fields, (Doctoral dissertation), The Joint Graduate School of Energy and Environment at King Mongkut's University of Technology Thonburi, Thailand. Smakgahn, K, Fumoto, T., & Yagi, K. (2009). Validation of revised DNDC model for methane emissions from irrigated rice fields in Thailand and sensitivity analysis of key factors. Journal of Geophysical Research, 114, 1-12. 272 Slocum, T. A., Mcmaster, R. B., Kessler, F. C., & Howard, H. H. (2009). Thematic Cartography and Geovisualization (3rd ed.).New Jersey, USA.: Prentice Hall Series in Geographic Information Science. Southeast Asia. (2009). In Encyclopædia Britannica. Retrieved April 5, 2009, from Encyclopedia Britannica Online: http://www.britannica.com/EBchecked/topic/556489/Southeast-Asia Sozanska, M., Skiba, U., & Metcalfe, S. (2002). Developing an inventory of N2O emissions from British soils. Atmospheric Environment, 36, 987– 998. Tan, B., Morisette, J.T., Wolfe, R.E., Gao, F., Ederer, G.A., Nightingale, J., & Pedelty, J.A. (2011). An enhanced TIMESAT algorithm for estimating vegetation phenology metrics from MODIS data. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. 4(2), 361-371. The University of Minnesota. (2011). MapServer open source web mapping. Retrieved from http://mapserver.org/ TMD (Thai Metorological Department). (2010). El Niño and La Niña Events. Research from Thai Metorological Department, Thailand. Tottrup, C., Rasmussen, M.S., Eklundh, L., & Jönsson, P. (2007). Mapping fractional forest cover across the highlands of mainland Southeast Asia using MODIS data and regression tree modeling. International Journal of Remote Sensing, 28(1), 23-46. Towprayoon, S., & Smakgahn, K. (2003). Study on model implementation on methane and nitrous oxide emissions from rice field. The Thailand Research Fund (TRF): research report. Towprayoon, S. (2005). Study of factors and their interaction on the effect of methane emission in rice field. The Thailand Research Fund (TRF): research report. Tsou, M.H. (2004). Integrated mobile GIS and wireless internet map servers for environmental monitoring and management. Cartography and Geographic Information Science, 31(3), 153-165. Tucker, C. J., Slayback, D. A., Pinzon, J. E., Los, S. O., Myneni, R. B., & Taylor, M. G. (2001). Higher northern latitude normalized difference vegetation index and growing season trends from 1982 to 1999. International Journal of Biometeorology, 45(4), 184-190. 273 United Nations Framework Convention on Climate Change (UNFCCC). (2007). Climate Change: Impacts Vulnerabilities and Adaptation in Developing Countries. Germany: UNFCCC. USDA. (2010). Thailand Grain and Feed Update Rice Update. Gain Report: Global Agricultural Information Network. Retrieved from http://gain.fas.usda.gov/Recent%20GAIN%20Publications/Grain%20and%20Feed%20U pdate_Bangkok_Thailand_12-23-2010.pdf Vietnam Trade Promotion Agency, Export Promotion Center. (2008). Report on Vietnamese Rice Sector. Retrieved from http://www.aseankorea.org/files/upload/board/120/6/VIETNAMESE%20RICE%20SECT OR%202008.pdf Wang, P.C.C. (1998). Graphical Representation of Multivariate Data. New York: Academic Press. Wardlow, B.D. Egbert, S.L., & Kastens, J.H. (2007). Analysis of time-series MODIS 250 m vegetation index data for crop classification in the U.S. Central Great Plains. Remote Sensing of Environment, 108, 290-310. Wassmann, R., & Dobermann, A. (2006). The 2nd Joint International Conference on sustainable energy and environment (SEE 2006), Bangkok, Thailand, 21-23 November 2006. White, M. A., Thornton, P. E., & Running, S.W. (1997). A Continental phenology model for monitoring vegetation responses to interannual climatic variability. Global Biogeochemical Cycles, 11(2), 217-234. White, M. A., & Nemani, R.R. (2006). Real-time monitoring and short-term forecasting of land surface phenology. Remote Sensing Of Environment, 104, 43-49. Xiao, X. M., Hagen, S., Zhang, Q. Y., Keller, M., & Moore, B. (2006a). Detecting leaf phenology of seasonally moist tropical forests in South America with multi-temporal MODIS images. Remote Sensing of Environment, 103(4), 465-473. Xiao, X. M., Boles, S., Frolking, S., Li, C. S., Babu, J. Y., Salas, W., et al. (2006b). Mapping paddy rice agriculture in South and Southeast Asia using multi-temporal MODIS images. Remote Sensing Of Environment, 100(1), 95-113. Xiao, X., Zhang, J., Yan, H., Wu, W., & Biradar, C. (2009). Land surface phenology: Convergence of satellite and CO2 eddy flux observation. In Noormets, A. (Ed.), Phenology of ecosystem processes (pp. 247-270). New York: Springer. 274 Xu, Z. X., Li,J. Y., & Liu, C. M. (2007). Long-term trend analysis for major climate variables in the Yellow River basin. Hydrological Processes, 21 (14), 1935-1948. Xu, H. (2011). A fine - scale analysis of nitrogen use efficency of paddy rice systems: a case study in Anji County in Zhejiang, China, (MSc thesis), Michigan State University. Young, K.B., Wailes, E.J., Cramer, G.L., & Khiem, N.T. (2002). Vietnam’s Rice Economy: Developments and Prospects. Retrieved from http://arkansasagnews.uark.edu/968.pdf Yu, F., Price, K.P., Ellis, J, & Shi, P. (2003). Response of seasonal vegetation development to climatic variations in eastern central Asia. Remote Sensing of Environment, 87, 42-54. Zhang, F, Li,C., Wang, Z., & Wu, H. (2006). Modeling impacts of management alternatives on soil carbon storage of farmland in Northwest China. Biogeosciences, 3, 451-466. Yue, S., Pilon, P., & Cavadias, G. (2002). Power of the Mann-Kendall and Spearman’s rho test for detecting monotonic trends in hydrological series. Journal of Hydrology, 259, 254271. Yue, H., Pilon, P., & Phinney, B. (2003). Canadian streamflow trend detection: impacts of serial and cross-correlation, Hydrological Sciences Journal, 48(1), 51-63 Zhang, X. Y., Hodges, J.C.F., Schaaf. C.B., Friedl, M. A., Strahler, A. H., & Gao, F. (2001).Global vegetation phenology from AVHRR and MODIS data. In Proceedings of the International Geoscience and Remote Sensing Symposium, University of New South Wales (Sydney, Australia: IGRSS) (on CD-ROM). Zhang, Y., Li, C., Zhou, X., & Moore, B. (2002). A simulation model linking crop growth and soil biogeochemistry for sustainable agriculture. Ecological Modelling, 151, 75-108. Zhang, X. Y., Friedl, M. A., Schaaf, C. B., Strahler, A. H., Hodges, J. C. F., Gao, F., et al. (2003). Monitoring vegetation phenology using MODIS. Remote Sensing of Environment, 84(3), 471-475. Zhang, X. Y., Friedl, M. A., Schaaf, C. B., & Strahler, A.H. (2005). Monitoring the response of vegetation phenology to precipitation in Africa by coupling MODIS and TRMM instruments. Journal of Geophysical Research,110(D12103), 1-14. Zhang, X.Y., Friedl, M. A., & Schaaf, C. B.. (2006). Global vegetation phenology from Moderate Resolution Imaging Spactroradiometer (MODIS): Evaluation of global patterns and comparison with in situ measurements. Journal of Geophysical Research,111(G04017), 1-14. 275 Zhang, L., Yu, D., Shi, X., Weindorf, D., Zhao, L., Ding, W., Wang, H., Pan, J., & Li, C. (2009). Quantifying methane emissions from rice fields in the Taihu Lake region, China by coupling a detailed soil database with biogeochemical model, Biogeosciences, 6, 739749. Zhang, Y., Wang, Y. Y., Su, S. L., & Li, C. S. (2011). Quantifying methane emissions from rice paddies in Northeast China by integrating remote sensing mapping with a biogeochemical model, Biogeosciences, 8, 1225-1235. Zhang, Y., Su, S., Zhang, F., Shi, R.,& Gao, W. (2012) Characterizing Spatiotemporal Dynamics of Methane Emissions from Rice Paddies in Northeast China from 1990 to 2010. Plos One, 7(1), 1-12. 276