Satellite-Based Estimation of Carbon Dioxide Budget in Tropical Peatland Ecosystems

: Tropical peatland ecosystems are known as large carbon (C) reservoirs and a ﬀ ect spatial and temporal patterns in C sinks and sources at large scales in response to climate anomalies. In this study, we developed a satellite data-based model to estimate the net biosphere exchange (NBE) in Indonesia and Malaysia by accounting for ﬁre emissions (FE), ecosystem respiration ( R e ), and gross primary production (GPP). All input variables originated from satellite-based datasets, e.g., the precipitation of global satellite mapping of precipitation (GSMaP), the land surface temperature (LST) of the moderate resolution imaging spectroradiometer (MODIS), the photosynthetically active radiation of MODIS, and the burned area of MODIS ﬁre products. First, we estimated the groundwater table ( GWT ) by incorporating LST and precipitation into the Keetch–Byram Drought Index (KBDI). The GWT was validated using in-situ measurements, with a root mean square error (RMSE) of 24.97 cm and an r-squared ( R 2 ) of 0.61. The daily GWT variations from 2002 to 2018 were used to estimate respiration ( R e ) based on a relationship between the in situ GWT and ﬂux-tower-observed R e . Fire emissions are a large direct source of CO 2 from terrestrial ecosystems into the atmosphere and were estimated by using MODIS ﬁre products and estimated biomass. The GPP was calculated based on the MODIS GPP product after parameter calibration at site scales. As a result, averages of long-term (17 years) R e , GPP, FE, and NBE from whole peatlands in the study area (6 ◦ N–11 ◦ S, 95–141 ◦ E) were 66.71, 39.15, 1.9, and 29.46 Mt C / month, respectively. We found that the NBE from tropical peatlands in the study area was greater than zero, acting as a C source. R e and FE were inﬂuenced by El Niño, and the value of the NBE was also high in the El Niño period. In future studies, the status of peatland degradation should be clariﬁed in detail to accurately estimate the C budget by applying appropriate algorithms of R e with delineations of types of anthropogenic impacts (e.g., drainages and ﬁres).


Introduction
Carbon dioxide (CO 2 ) is one of the gases causing global warming [1]. Peatland ecosystems are expected to play an essential role in global warming mitigation as carbon reservoirs [2]. The impacts of peatlands on greenhouse gas emissions have been examined [3,4]. Despite covering only 3% of the global land area, peatlands contain about one-third of global soil organic carbon [5]. Peat soils contain a high carbon density and consist of undecomposed carbon that is induced by a large amount of biomass under shallow groundwater table (GWT) conditions [6]. Recent disturbances in the ecosystem due to deforestation, wildfires, and groundwater drainages have significantly altered their CO 2 budget [3]. Therefore, estimating their carbon budget and proposing methods to suppress CO 2 emission across from MODIS; (3) to calibrate the satellite-based GPP of MODIS by fitting it to the flux tower-observed GPP; and (4) to integrate emission (Re, FE) components and an absorption (GPP) component into a CO2 budget of target peatland areas.

Study Area
Southeast Asia (SE Asia) region is known as a hotspot for tropical peatlands, accounting for half of the global peatland emissions due to rapid deforestation and drainage [14]. Within SE Asia, Indonesia has the highest CO2 emission rate due to peatland decomposition and fires. The target area for this study was the peatlands in Indonesia and Malaysia in the Borneo Island region (6°N-11°S, 95-141°E) (Figure 1). In the study area, the total peat forest area was 205,408 km 2 , as reported by the Global Forest Watch collected from the Ministry of Agriculture (2012) and Wetlands International (2004) for Indonesia and Malaysia [35,36]. The locations of flux towers are shown in Figure 1. The information of latitude and longitude is as follows: undrained forest (UF) (2°20′42″S, 114°2′11″E), drained forest (DF) (2°20′59″S, 114°8′13″E), and burned forest (BF) (2°20′14″S, 114°2′13″E).  [36] was used as a peatland map. The locations of the flux towers are marked as red circles. UF, DF, and BF denote undrained forest, drained forest, and burned forest, respectively.

KBDI-Based GWT
The GWT was an essential factor for estimating Re in this study. To estimate the spatially distributed GWT, we used satellite-based data in the computation of the KBDI. As the input data, we used LST from a multi-functional transport satellite (MTSAT) (4 km spatial resolution) and precipitation (P) from global satellite mapping of precipitation (GSMaP) (10 km spatial resolution). Since no MTSAT-based LST data in 2002-2006 were available due to the satellite's pre-observation period, we instead used MODIS LST (MOD11A1, 1 km spatial resolution) for these years. MTSAT is a geostationary satellite. MODIS sensors are onboard a polar-orbiting Terra satellite. Therefore, their observation frequencies are different. The maximum LST (LSTmax) was calculated by using the hourly LST by MTSAT data and the daytime product by Terra MODIS data. The merit of the use of MTSAT LST is that an accurate LSTmax could be obtained due to its high temporal observation. However, for analyzing the relationship between climatic factors and Re, the long-term interannual variations were also crucial. Therefore, the MODIS LST altered the MTSAT one. As an ancillary dataset, the land surface water coverage (LSWC) of the Advanced Microwave Scanning Radiometer for Earth Observation System (AMSR-E) was used with a 10 km spatial resolution for confirming the existence of surface water [21]. The LSWC was calculated by combining MODIS normalized Figure 1. Spatial distribution of tropical peatlands in this study (colored by brown). The global forest watch database (GFW) [36] was used as a peatland map. The locations of the flux towers are marked as red circles. UF, DF, and BF denote undrained forest, drained forest, and burned forest, respectively.

KBDI-Based GWT
The GWT was an essential factor for estimating R e in this study. To estimate the spatially distributed GWT, we used satellite-based data in the computation of the KBDI. As the input data, we used LST from a multi-functional transport satellite (MTSAT) (4 km spatial resolution) and precipitation (P) from global satellite mapping of precipitation (GSMaP) (10 km spatial resolution). Since no MTSAT-based LST data in 2002-2006 were available due to the satellite's pre-observation period, we instead used MODIS LST (MOD11A1, 1 km spatial resolution) for these years. MTSAT is a geostationary satellite. MODIS sensors are onboard a polar-orbiting Terra satellite. Therefore, their observation frequencies are different. The maximum LST (LST max ) was calculated by using the hourly LST by MTSAT data and the daytime product by Terra MODIS data. The merit of the use of MTSAT LST is that an accurate LST max could be obtained due to its high temporal observation. However, for analyzing the relationship between climatic factors and R e , the long-term interannual variations were also crucial. Therefore, the MODIS LST altered the MTSAT one. As an ancillary dataset, the land surface water coverage (LSWC) of the Advanced Microwave Scanning Radiometer for Earth Observation System (AMSR-E) was used Remote Sens. 2020, 12, 250 4 of 21 with a 10 km spatial resolution for confirming the existence of surface water [21]. The LSWC was calculated by combining MODIS normalized difference water index (NDWI) and the AMSR-E daily normalized difference polarization index (NDPI) [37].

Flux Tower Observations
The study area contained three sites: an undrained forest (UF), a drained forest (DF), and a drained burnt forest (BF) on flat terrain within 15 km [25] in Palangkaraya, Kalimantan, Indonesia (Palangkaraya drained forest (PDF) site in AsiaFlux database). The available observation periods are as follows: DF was from 2001, and UF and BF were from 2004. The heights of the flux towers are 35, 50, and 4 m for UF, DF, and BF, respectively. The DF and BF towers were drained and located about 400 and 200 m from the canals [25]. UF is a relatively intact forest site located in a national park along the Sebangau River. The compositions of land cover types in the UF, DF, and BF based on MODIS land cover types (MCD12Q1) by the international geosphere-biosphere programme (IGBP) scheme were evergreen broadleaf forest (EBF) (in UF and DF), woody savannas (WSAV) (in DF), and grasslands (GL) and savannas (SAV) (in BF). The spatial representativeness of flux tower is usually known as 100-1000 m [38]. CO 2 flux data were generated every 30 min. Each 5-min CO 2 concentration that was recorded with a closed-path CO 2 analyzer (LI7500, Li-Cor Inc., Lincoln, NE, USA) was interpolated for estimating the 30-min CO 2 fluxes. The observed GWT was also used and was observed in 30-min intervals. Flux partitioned R e and the observed GWT were used for modeling R e . The GWT was used for the validation of the RS-based GWT. The flux partitioned GPP was also used to calibrate parameters in the MODIS GPP algorithm.

MODIS GPP
The GPP is the largest carbon flux into terrestrial ecosystems. Many peatlands include forests in the area. The PDF flux tower site also has a large forested area of~50 km 2 . As usual for tropical forests, the photosynthesis in this area is high (more than 10 g C/m 2 /day). The 8-day MODIS GPP product (MOD17A2H) with 500 m resolution was used to estimate from the regional to national scales of GPPs for a long-term period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). MOD17A2H (Collection 6) has improved its spatial resolution by using high-resolution reanalysis meteorological data (modern era-retrospective-analysis for research and applications (MERRA)-2 of Global Modeling and Assimilation Office of the National Aeronautics and Space Administration (NASA)) [39]. The MODIS GPP algorithm is a typical biochemical model that uses light use efficiency for each biome type based on the land cover map (MCD12Q1). Temperature and water stresses are the factors suppressing the potential GPP for the vegetation type. Absorbed photosynthetically active radiation (APAR) was calculated with solar radiations of reanalysis data, and the fraction of APAR (fAPAR) was multiplied with the incident PAR (IPAR) on the vegetative surface. Uncertainties in the data exist due to low PAR under the cloud-covered situations [40]. The MODIS sensor was installed on a Terra satellite, and its revisit time was every 10:30 in local time [34].

Fire Emissions Based on MODIS
Satellite-based thermal information is suitable to detect fire events. We estimated CO 2 emissions from anthropogenic-induced and wildfires as another flux in the calculation of CO 2 balance. During a fire event, various kinds of trace gases are released [41]. In this study, we focused on CO 2 emissions. CO 2 emissions from fire events were estimated by using biomass, burned areas, fuel types, combustion severities, and emission factors. For biomass, we used the vegetation integrative simulator for trace gases (VISIT), a process-based model, with 0.5 • spatial and monthly temporal resolutions [28,42,43]. In the VISIT, aboveground biomass (AGB) and SOM are calculated through process-based algorithms of the terrestrial carbon cycle. The MODIS normalized difference vegetation index (NDVI) (a 1 km spatial resolution) was used to downscale the VISIT-based biomass estimation into a 1 km resolution by using logistic regression functions for each IGBP biome type in MCD12Q1 (Table A1) [44]. Burned areas and their dates were obtained from a MODIS fire product, MCD45A1 [45]. Fuel types were Remote Sens. 2020, 12, 250 5 of 21 obtained from the MCD12Q1 500 m land cover map. Differences in the NDVI before and after fire events were directly reflected the loss of AGB. Combustion severities were represented by fire radiative power (FRP) in MODIS fire products. Wooster et al. [46] derived a combustion completeness ratio by using FRP. The merits of using FRP include its ability to describe daily temporal dynamics of fire severity. In addition, its spatial resolution is high (1 km). Emission factors were obtained from the global fire emissions database (GFED) 4.1s [47], meaning the proportion of CO 2 emissions within biome-specific carbon gases by biomass burning. For the CO 2 emission ratios of AGB and SOM, 0.49 and 0.57 were applied, respectively [33].

GWT Calculations
The KBDI was originally developed by Keetch and Byram in 1968 for a wildfire monitoring system [20]. It explains how the interaction between P (mm) and air temperatures (T air , • C) changes and how the evapotranspiration differs in specific vegetation types interpreting in situ measurements from various fire events: where dQ is a drought factor, Q 0 is the initial moisture deficiency (mm) at time zero, T air is air temperature, P is precipitations, and dτ is a time interval (1 day in this study). The index ranges from 0 (moist) to 800 (dry) (unitless) [20]. The range of the index is determined by assuming that fully saturated soil has 8 inches of moisture readily available to vegetation [21]. The original index was derived from field observations. For further application to satellite-based data, a previous study-modified KBDI referring to in situ measurements of the GWT [21] and the equations were used in this study: where mKBDI is the modified KBDI for using satellite-based data; mKBDI0 is an initial value that was calibrated with in situ measurements of the GWT, p day is a daily precipitation (mm/day), LST max is maximum daily LST (degrees in Celsius), P ann is an annual precipitation in each year (mm/year), and the previous study [21] added it to account for spatial variations in regional climatic moisture conditions. The initial value (mKBDI0) was determined by comparing a consistency of estimated inundation timing (i.e., GWT > 0) and the AMSR-E inundation map. CO 2 emissions are associated with changes in local hydrology, such as fluctuations in groundwater levels. CO 2 emissions from peatlands are highly related to the GWT [25], and this study used the GWT as an explanatory variable of R e ( Figure 2). Exponentially regressed p day , LST max , and P ann are factors that represent the daily water recharge, evaporation, and climatic long-term water recharge, respectively. In this study, the T air of original KBDI (Equation (1)) was altered by the satellite-based LST in the mKBDI (Equation (2)), since it is straightforward to use satellite-based LST, and the usefulness has been evaluated for substituting with T air , showing high correlation coefficient over 0.80 (r 2 ) [48,49]. A linear regression function, established by the relationship between the satellite-based daily mKBDI and in situ GWT measurements, was used to estimate the GWT spatial variation (Equation (3)) [21]. The AMSR-E-based LSWC, which represents land surface inundation, was used as supplemental data to generate the GWT map [21]. If an LSWC existed and the modelled GWT showed above 0 m, the initial value of the GWT of the day was set to 0.
Remote Sens. 2019, 11, x FOR PEER REVIEW 2 of 3 surface moisture condition was relatively rapid and simple once we considered that the KBDI-GWT represents a shallow GWT [50].

Ecosystem Respiration
Re was calculated by using Equations (4)- (6). These equations were determined from observational data in a previous study [25] and the data distribution was visualized ( Figure 3). Flux tower observations measure the total carbon budget at heights more than 40 m above the land surface. The eddy covariance method is used to determine Re and the GPP by basically calculating daytime and nighttime carbon fluxes. In the flux tower sites of this study, observed half-hourly nighttime net ecosystem CO2 exchange (NEE) was used as nighttime Re. A look-up-table (LUT) method combining environmental factors (e.g., soil moisture below 5 cm and GWL) was applied to model nighttime Re. For the next step, nighttime Re was extrapolated to daytime Re by using the established model. The estimated Re was used to infer the daytime GPP by combining an LUT with PAR [25]. Finally, Re in this study represents the CO2 emissions of the biological activity of plants including soil respiration [25,51]. By using Equations (4)-(6) for UF, DF, and BF, the Re values from whole tropical peatlands in Indonesia and Borneo of Malaysia were simulated. Furthermore, we discussed the potential impacts of peatlands conditions that had undergone drainage and fires. Finally, we used the equation of Re (UF) (Equation (4)) when the Re values of the whole peatland area were calculated, assuming that all peatlands from the GFW (reference map) were intact (natural) forest (e.g., evergreen broadleaf forest (for more land cover types information of DF and BF, refer to Section 2.2.2)). Notably, we assumed that all peatland regions in this study were dense forest since a relevant dataset for describing the spatial distribution of disturbance levels (drainages, fires, and post-fire regrowth) was not available. Flowchart of this study. As CO 2 fluxes emit to the atmosphere, ecosystem respiration (R e ), and fire emissions (FE) were calculated with groundwater table (GWT) and moderate resolution imaging spectroradiometer (MODIS) fire products, respectively. As a CO 2 flux from the atmosphere into plants, the gross primary production (GPP) of the MODIS product was used. Especially, the MODIS GPP was calibrated via flux tower observation in Palangkaraya in Indonesia. Finally, the net biosphere exchange (NBE) was calculated by balancing R e , FE, and GPP. Dashed rectangles indicate corresponding sections in this paper.
A bilinear method was used to generate the 1 km GWT data from P (10 km) and LST (4 km). The discrepancy of spatial resolutions between the input data sources might cause the presence of mixed information at the sub-pixel scale in complex terrain, e.g., mountainous areas. Notably, most of the target areas in this study are located coastal or lowland regions, and the response of the GWT to the surface moisture condition was relatively rapid and simple once we considered that the KBDI-GWT represents a shallow GWT [50].

Ecosystem Respiration
R e was calculated by using Equations (4)-(6). These equations were determined from observational data in a previous study [25] and the data distribution was visualized ( Figure 3). Flux tower observations measure the total carbon budget at heights more than 40 m above the land surface. The eddy covariance method is used to determine R e and the GPP by basically calculating daytime and nighttime carbon fluxes. In the flux tower sites of this study, observed half-hourly nighttime net ecosystem CO 2 exchange (NEE) was used as nighttime R e . A look-up-table (LUT) method combining environmental factors (e.g., soil moisture below 5 cm and GWL) was applied to model nighttime R e . For the next step, nighttime R e was extrapolated to daytime R e by using the established model. The estimated R e was used to infer the daytime GPP by combining an LUT with PAR [25]. Finally, R e in this study represents the CO 2 emissions of the biological activity of plants including soil respiration [25,51]. By using Equations (4)-(6) for UF, DF, and BF, the R e values from whole tropical peatlands in Indonesia and Borneo of Malaysia were simulated. Furthermore, we discussed the potential impacts of peatlands conditions that had undergone drainage and fires. Finally, we used the equation of R e (UF) (Equation (4)) when the R e values of the whole peatland area were calculated, assuming that all peatlands from the GFW (reference map) were intact (natural) forest (e.g., evergreen broadleaf forest (for more land cover types information of DF and BF, refer to Section 2.2.2)). Notably, we assumed that all peatland regions in this study were dense forest since a relevant dataset for describing the spatial distribution of disturbance levels (drainages, fires, and post-fire regrowth) was not available.

GPP Calibration
The MODIS-based GPP (MOD17A2H) time series underestimated the GPP derived from flux towers when it was validated with the GPP of flux towers for three conditions (i.e., UF, DF, and BF) in Palangkaraya. In all cases, the MODIS GPP underestimated flux tower one, and the R 2 in DF was around 0.29 based on the least squares regression [52] between the MODIS GPP and the flux towerderived GPP (Equation (7), and Figure 4).
It has also been reported that the least squares regression method is sensitive to outliers [52]. To overcome this problem, we applied the robust regression method. This method has the following advantages. The robust regression method judges whether a point is an outlier or not. The coefficients of regression function are calculated by minimizing the absolute error. Furthermore, a weight function derived by the data distribution was applied to remove the excessive effects that were induced by outliers [53]. Hence, the weight function of each data point in a robust regression is useful to represent a best fitting line. Consequently, the linear equation with the robust regression method can reflect a more robust data distribution, removing the effects of outliers. In this study, some ranges in the MOD17A2H GPP, when the flux tower GPP was from 4 to 6 g C/m 2 /day, distributed likely outliers on the scatter plot ( Figure 4). Underestimation in the low-end ranges of the MOD17A2H GPP compared with the flux tower GPP (e.g., the MOD17A2H GPP showed almost 0 g C/m 2 /day, but the flux tower GPP was more than 4 g C/m 2 /day), might have been caused by diffuse PAR [54]. Diffuse PAR generally provides more photosynthesis activity than that of direct PAR, and the MOD17 GPP algorithm does not account for this effect. On the other hand, higher values in the MODIS than flux the tower GPP could be produced by an overestimation of the MOD17 GPP algorithm. A clear sky condition with excess solar radiation might increase the potential GPP of MODIS, which is induced by incident PAR; however, the real photosynthesis activity would be saturated [55]. If a robust regression method [53] was applied to generate the linear regression for the calibration of the MODIS GPP, the R 2 would be increased around 0.37 (Equation (8), and Figure 4). A function in MATLAB software named the "robustdemo" was used to perform a robust regression method.
The two regressions were applied to the empirical parameter calibration in this study. As a result, the linear fitting function of robust regression (Figure 4) seemed to explain more plots than

GPP Calibration
The MODIS-based GPP (MOD17A2H) time series underestimated the GPP derived from flux towers when it was validated with the GPP of flux towers for three conditions (i.e., UF, DF, and BF) in Palangkaraya. In all cases, the MODIS GPP underestimated flux tower one, and the R 2 in DF was around 0.29 based on the least squares regression [52] between the MODIS GPP and the flux tower-derived GPP (Equation (7), and Figure 4).
It has also been reported that the least squares regression method is sensitive to outliers [52]. To overcome this problem, we applied the robust regression method. This method has the following advantages. The robust regression method judges whether a point is an outlier or not. The coefficients of regression function are calculated by minimizing the absolute error. Furthermore, a weight function derived by the data distribution was applied to remove the excessive effects that were induced by outliers [53]. Hence, the weight function of each data point in a robust regression is useful to represent a best fitting line. Consequently, the linear equation with the robust regression method can reflect a more robust data distribution, removing the effects of outliers. In this study, some ranges in the MOD17A2H GPP, when the flux tower GPP was from 4 to 6 g C/m 2 /day, distributed likely outliers on the scatter plot ( Figure 4). Underestimation in the low-end ranges of the MOD17A2H GPP compared with the flux tower GPP (e.g., the MOD17A2H GPP showed almost 0 g C/m 2 /day, but the flux tower GPP was more than 4 g C/m 2 /day), might have been caused by diffuse PAR [54]. Diffuse PAR generally provides more photosynthesis activity than that of direct PAR, and the MOD17 GPP algorithm does not account for this effect. On the other hand, higher values in the MODIS than flux the tower GPP could be produced by an overestimation of the MOD17 GPP algorithm. A clear sky condition with excess solar radiation might increase the potential GPP of MODIS, which is induced by incident PAR; however, the real photosynthesis activity would be saturated [55]. If a robust regression method [53] was applied to generate the linear regression for the calibration of the MODIS GPP, the R 2 would Remote Sens. 2020, 12, 250 8 of 21 be increased around 0.37 (Equation (8), and Figure 4). A function in MATLAB software named the "robustdemo" was used to perform a robust regression method. robust regression method (Equation (8) where Modified GPPleast and Modified GPProbust are the results from modification with the flux tower GPP by using the least squares and the robust regression method, respectively, and MOD17A2H GPP is results of the MODIS GPP product ( Figure 4). Finally, the Modified GPProbust (Equation (8)) was used for estimating the NBE in whole tropical peatland area in this study. We used a model established at the DF flux tower because the GPP in DF showed the highest R 2 among the three towers.

Net Biosphere Exchange
The NBE is determined from Re and FE as carbon sources and the GPP as a carbon absorption (Equation (9), Figure 2). To effectively conduct validation within mixtures of positive (+) and negative (-) values of CO2 fluxes, the RMSE of the NBE was normalized (NRMSE) by using ymax -ymin (Equation (10)) where y is the result of Re, FE, GPP, or NBE in this study. The two regressions were applied to the empirical parameter calibration in this study. As a result, the linear fitting function of robust regression (Figure 4) seemed to explain more plots than that of the least squares regression method. The root mean squared error (RMSE) and standard deviation values of the original MODIS GPP data, Modified GPP least , and Modified GPP robust were 3.91 ± 2.95, 1.47 ± 1.59, and 1.33 ± 1.40 g C/m 2 /day (± standard deviation), respectively. We concluded the linear fitting equation of robust regression was more suitable for explaining flux tower GPP when using MOD17A2H. Based on the comparison, we modified the MODIS GPP by using the result of the robust regression method (Equation (8));

GWT
Modi f ied GPP robust = 2.1142(MOD17A2H GPP) − 12.589 (8) where Modified GPP least and Modified GPP robust are the results from modification with the flux tower GPP by using the least squares and the robust regression method, respectively, and MOD17A2H GPP is results of the MODIS GPP product (Figure 4). Finally, the Modified GPP robust (Equation (8)) was used for estimating the NBE in whole tropical peatland area in this study. We used a model established at the DF flux tower because the GPP in DF showed the highest R 2 among the three towers.

Net Biosphere Exchange
The NBE is determined from R e and FE as carbon sources and the GPP as a carbon absorption (Equation (9), Figure 2). To effectively conduct validation within mixtures of positive (+) and negative Remote Sens. 2020, 12, 250 9 of 21 (-) values of CO 2 fluxes, the RMSE of the NBE was normalized (NRMSE) by using y max -y min (Equation (10)) where y is the result of R e , FE, GPP, or NBE in this study.

GWT
Daily GWT averages from the KBDI and the in-situ measurements in UF, DF, and BF were 23.5, 14.72, 55.6, and 12.18 cm, respectively. RMSEs between the KBDI-based and in-situ GWTs were 24.97, 39.32, and 21.03 cm, respectively. The largest discrepancy was shown in DF among the three sites, since the GWT was artificially decreased through drainage canals. When only UF and BF were considered, the RMSE of the KBDI-based and in-situ GWTs was 20.65 cm. Inundated seasons, which showed a GWT greater than zero, usually appeared from November to April. As results of linear fitting between the KBDI-based and in-situ GWTs, the R 2 values of the UF, DF, and BF sites were 0.61, 0.64, and 0.66, respectively (Figure 5b-d). In terms of the slope of linear fittings, the GWTs at the UF site were closer to a 1:1 relationship. In contrast, the in-situ GWT at the BF site showed less variations in data distributions than the variations of KBDI-based GWT with a slope of the linear function of 0.77. The R 2 of BF was the highest; both the GWTs matched well under very dry (low GWT) conditions. Since BF had undergone disturbance events (e.g., wildfires) and less vegetation cover than other sites, the characteristics of evapotranspiration could be simply explained. Daily GWT averages from the KBDI and the in-situ measurements in UF, DF, and BF were 23.5, 14.72, 55.6, and 12.18 cm, respectively. RMSEs between the KBDI-based and in-situ GWTs were 24.97, 39.32, and 21.03 cm, respectively. The largest discrepancy was shown in DF among the three sites, since the GWT was artificially decreased through drainage canals. When only UF and BF were considered, the RMSE of the KBDI-based and in-situ GWTs was 20.65 cm. Inundated seasons, which showed a GWT greater than zero, usually appeared from November to April. As results of linear fitting between the KBDI-based and in-situ GWTs, the R 2 values of the UF, DF, and BF sites were 0.61, 0.64, and 0.66, respectively (Figure 5b-d). In terms of the slope of linear fittings, the GWTs at the UF site were closer to a 1:1 relationship. In contrast, the in-situ GWT at the BF site showed less variations in data distributions than the variations of KBDI-based GWT with a slope of the linear function of 0.77. The R 2 of BF was the highest; both the GWTs matched well under very dry (low GWT) conditions. Since BF had undergone disturbance events (e.g., wildfires) and less vegetation cover than other sites, the characteristics of evapotranspiration could be simply explained.

Ecosystem Respiration
The RMSE and bias of satellite-based Re were 1.12 and 0.09 g C/m 2 /day, respectively, for the all three sites. In detail, the RMSEs of UF, DF, and BF were 1.31, 1.22, and 0.82 g C/m 2 /day, respectively. The biases were −0.15, 0.45, and −0.02 g C/m 2 /day for UF, DF, and BF, respectively. The standard deviations were 1.14, 0.8, and 0.5, respectively. As shown in Figure 6, the variances of flux tower-

Ecosystem Respiration
The RMSE and bias of satellite-based R e were 1.12 and 0.09 g C/m 2 /day, respectively, for the all three sites. In detail, the RMSEs of UF, DF, and BF were 1.31, 1.22, and 0.82 g C/m 2 /day, respectively. The biases were −0.15, 0.45, and −0.02 g C/m 2 /day for UF, DF, and BF, respectively. The standard deviations were 1.14, 0.8, and 0.5, respectively. As shown in Figure 6, the variances of flux tower-based R e were highest in UF and lowest in BF, which was related to their vegetation coverages. UF and DF are relatively dense forests; however, BF has lost its vegetation due to fires, resulting in a sparse vegetation area. DF is a reforested area that is covered by secondary forests for logging [25]. The GPP values of DF and BF were obviously lower than those of the intact site (UF). Natural forests usually have more dense vegetation coverages. Vegetation productivities are related to the amount of biomass, including leaves, foliage, branches, plant body, and litter. Consequently, R e is also high in natural forests. The RMSE, bias and standard deviation of UF were higher than the others. The NRMSEs for UF, DF, and BF were 0.01, 0.01, and 0.02, respectively. From this result, BF showed more uncertainties in Re, when considering original data distributions, than the other sites. For DF, a trend of a lower GWT caused by drainage was already confirmed in Section 3.1. The lower GWT of DF than that of UF had the effect of increasing Re, as well. The in situ Re of DF showed a high positive bias (0.45 g C/m 2 /day) due to drier conditions that promoted the decomposition of peat soil. Note the KBDI-based GWT for Re calculation in DF did not consider the effects of drainage canals. Hence, we concluded that UF provided the most suitable model to extrapolate Re flux observations to spatial estimations. From the annual average maps of Re from 2002 to 2018, the carbon (C) emissions by Re were high (814.63 Mt C/year), middle (754.76 Mt C/year), and low (404.70 Mt C/year), in the UF, DF, and BF sites, respectively (Figure 7). Pixels with values larger than 1000 g C/m 2 /year were not confirmed in the Re maps. This means that the method proposed in this study to estimate Re only depends on a formula that represents the groundwater condition in peatlands. Choosing appropriate ground conditions is necessary and affects spatial distribution of Re. However, we only considered a baseline of Re when using the data at the UF site, assuming that all peatlands were intact and the Re was high due to the forest productivity and autotrophic respiration being high, since those are related to vegetation growth rate [56].  The RMSE, bias and standard deviation of UF were higher than the others. The NRMSEs for UF, DF, and BF were 0.01, 0.01, and 0.02, respectively. From this result, BF showed more uncertainties in R e , when considering original data distributions, than the other sites. For DF, a trend of a lower GWT caused by drainage was already confirmed in Section 3.1. The lower GWT of DF than that of UF had the effect of increasing R e , as well. The in situ R e of DF showed a high positive bias (0.45 g C/m 2 /day) due to drier conditions that promoted the decomposition of peat soil. Note the KBDI-based GWT for R e calculation in DF did not consider the effects of drainage canals. Hence, we concluded that UF provided the most suitable model to extrapolate R e flux observations to spatial estimations. From the annual average maps of R e from 2002 to 2018, the carbon (C) emissions by R e were high (814.63 Mt C/year), middle (754.76 Mt C/year), and low (404.70 Mt C/year), in the UF, DF, and BF sites, respectively (Figure 7). Pixels with values larger than 1000 g C/m 2 /year were not confirmed in the R e maps. This means that the method proposed in this study to estimate R e only depends on a formula that represents the groundwater condition in peatlands. Choosing appropriate ground conditions is necessary and affects spatial distribution of R e . However, we only considered a baseline of R e when using the data at the UF site, assuming that all peatlands were intact and the R e was high due to the forest productivity and autotrophic respiration being high, since those are related to vegetation growth rate [56].

GPP in Tropical Peatlands
The daily average and the standard deviation of the modified GPProbust were 8.91 and 1.

GPP in Tropical Peatlands
The daily average and the standard deviation of the modified GPP robust were 8.91 and 1.21 g C/m 2 /day, respectively, at the Palangkaraya DF site from 2002 to 2018 (Figure 8). The original data from the DF flux tower were 9.5 and 0.74 g C/m 2 /day for the daily average and the standard deviation, respectively. Before the modification of the GPP in this study, although the daily average was low (6.26 g C/m 2 /day), the standard deviation of the MODIS GPP was relatively high (2.57 g C/m 2 /day). The area's average of the monthly GPP was 39.15 Mt C/month, and its standard deviation was 1.09 Mt C/month). The mean annual GPP ranged from 2100 to 2400 g C/m 2 /year (Figure 9).

Fire Emissions
The monthly average FE from 2002 to 2018 was 1.90 Mt C/month in the study area. The FE from SOM that represented the C emissions from burning soil organic biomass was 1.6 Mt C/month on average, which accounted for 84.17% of the total FE including AGB and SOM. FE from AGB was 0.3 Mt C/month and accounted for 15.83% of the total FE. In this study area, FE from SOM was the main contributor to the total FE. Monthly time series of FE and GFED were compared (Figure 10a). The R 2 between both datasets was 0.74, the RMSE was 1.67 Mt C/month, and the slope of linear fitting function was 0.80 (Figure 10b). The spatial distributions of the average of yearly FE are shown in Figure 11. As a result, the FE values were higher in Sumatra and Kalimantan of Indonesia and Borneo of Malaysia than the other regions.

Fire Emissions
The monthly average FE from 2002 to 2018 was 1.90 Mt C/month in the study area. The FE from SOM that represented the C emissions from burning soil organic biomass was 1.6 Mt C/month on average, which accounted for 84.17% of the total FE including AGB and SOM. FE from AGB was 0.3 Mt C/month and accounted for 15.83% of the total FE. In this study area, FE from SOM was the main contributor to the total FE. Monthly time series of FE and GFED were compared (Figure 10a). The R 2 between both datasets was 0.74, the RMSE was 1.67 Mt C/month, and the slope of linear fitting function was 0.80 (Figure 10b). The spatial distributions of the average of yearly FE are shown in Figure 11. As a result, the FE values were higher in Sumatra and Kalimantan of Indonesia and Borneo of Malaysia than the other regions.

Fire Emissions
The monthly average FE from 2002 to 2018 was 1.90 Mt C/month in the study area. The FE from SOM that represented the C emissions from burning soil organic biomass was 1.6 Mt C/month on average, which accounted for 84.17% of the total FE including AGB and SOM. FE from AGB was 0.3 Mt C/month and accounted for 15.83% of the total FE. In this study area, FE from SOM was the main contributor to the total FE. Monthly time series of FE and GFED were compared (Figure 10a). The R 2 between both datasets was 0.74, the RMSE was 1.67 Mt C/month, and the slope of linear fitting function was 0.80 (Figure 10b). The spatial distributions of the average of yearly FE are shown in

NBE Results and Validation with Flux Tower Data
The NBE of tropical peatlands was calculated based on Equation (9). By using the Re model established at the UF site, the monthly average and the standard deviation of the NBE were estimated as 29.46 and 6.26 Mt C/month, respectively. When the Re values of BF and DF were applied, the NBE values were estimated as −4.57 ± 4.12 and 24.91 ± 5.04 Mt C/month on average (± standard deviations), respectively. From this estimation, the monthly average GPP (39.15 Mt C/month) was equivalent to 58.68% of Re (66.71 Mt C/month), and the FE (1.90 Mt C/month) accounted for 2.84% of Re. Consequently, the NBE and resultant emissions of the monthly time series (Figure 12) were always high.

NBE Results and Validation with Flux Tower Data
The NBE of tropical peatlands was calculated based on Equation (9). By using the Re model established at the UF site, the monthly average and the standard deviation of the NBE were estimated as 29.46 and 6.26 Mt C/month, respectively. When the Re values of BF and DF were applied, the NBE values were estimated as −4.57 ± 4.12 and 24.91 ± 5.04 Mt C/month on average (± standard deviations), respectively. From this estimation, the monthly average GPP (39.15 Mt C/month) was equivalent to 58.68% of Re (66.71 Mt C/month), and the FE (1.90 Mt C/month) accounted for 2.84% of Re. Consequently, the NBE and resultant emissions of the monthly time series (Figure 12) were always high.

NBE Results and Validation with Flux Tower Data
The NBE of tropical peatlands was calculated based on Equation (9). By using the R e model established at the UF site, the monthly average and the standard deviation of the NBE were estimated as 29.46 and 6.26 Mt C/month, respectively. When the R e values of BF and DF were applied, the NBE values were estimated as −4.57 ± 4.12 and 24.91 ± 5.04 Mt C/month on average (± standard deviations), respectively. From this estimation, the monthly average GPP (39.15 Mt C/month) was equivalent to 58.68% of R e (66.71 Mt C/month), and the FE (1.90 Mt C/month) accounted for 2.84% of R e . Consequently, the NBE and resultant emissions of the monthly time series (Figure 12) were always high. as 29.46 and 6.26 Mt C/month, respectively. When the Re values of BF and DF were applied, the NBE values were estimated as −4.57 ± 4.12 and 24.91 ± 5.04 Mt C/month on average (± standard deviations), respectively. From this estimation, the monthly average GPP (39.15 Mt C/month) was equivalent to 58.68% of Re (66.71 Mt C/month), and the FE (1.90 Mt C/month) accounted for 2.84% of Re. Consequently, the NBE and resultant emissions of the monthly time series (Figure 12) were always high.

Spatial Distributions of GPP and Potential Uncertainties
We observed small regional variations in the GPP (Figure 9). Similarly to R e distributions, meteorological conditions across the peatlands in the study area were evenly distributed, since the input variables of the MODIS GPP algorithm relied on a minimum of daily temperature (TMIN), vapor pressure deficit (VPD), and fPAR [39,40]. Underestimations of the MOD17A2H GPP compared with the flux tower-based GPP were related to uncertainties of input variables in the GPP algorithm. In particular, MODIS fPAR data were essential because they were used to estimate detailed spatial and temporal patterns in the vegetation status. However, the data are acquired once per day at 10:30 am (local time). This may have caused a discrepancy between the MODIS GPP and the flux tower-based GPP. The other input variables of the MODIS algorithm were obtained from reanalysis datasets (e.g., MERRA-2) through numerical calculations of, e.g., solar radiation and relative humidity. Since reanalysis datasets have a coarser spatial resolution (~12.5 km) than the GPP product, resampling methods could have caused additional uncertainties [39,40]. The accurate estimation of near-surface relative humidity or VPD is also challenging [57,58]. Discrepancies between fPAR and the real surface situation including land cover changes could have occurred, since the land cover maps are updated on an annual basis (MCD12Q1); however, flux tower observations provide daily or shorter-interval GPP.
Some uncertainties could have been induced by the discrepancy between the spatial representativeness of the flux towers and the MOD17A2H GPP. As mentioned above, the footprints of flux towers vary (100-1000 m); however the MOD17A2H GPP shows a fixed spatial resolution (500 m). Moreover, changes in sub-pixel scales might not be detectable for the MOD17A2H GPP [59], although flux towers can capture them. In addition, the spatial heterogeneity in both biome types could have resulted in the discrepancy in the process of calibration for whole peatland areas. For future studies, with more flux towers representing peatlands and associated biomes, the GPP calibration could be improved by considering a light use efficiency (LUE) and plant functional types by using satellite-based land cover maps [60].
We also compared our results with another study, the FLUXNET multi-tree ensemble (MTE) based GPP [24]. As a result of the comparison, averages and standard deviations were 33.07 ± 39.11 and 1.21 ± 1.10 (Tg C/month), respectively, for FLUXNET-MTE and this study. The RMSE and bias were 6.31 and 6.03 (Tg C/month), respectively. Though the spatial resolutions of both GPP datasets differed by 100 times, the amount of the GPP and the seasonal patterns between both datasets were demonstrably similar (Figure is not shown).

FE Comparisons with GFED 4.1s
Both the FE of this study and the GFED 4.1s showed a good agreement on fluctuation patterns and their magnitude (Figure 10a). Higher GFED values than those of FE in general cases were potentially due to the difference in spatial resolution (GFED has a 25 km resolution and that of FE is 1 km). Another potential discrepancy was caused by the calculation of combustion completeness. GFED uses fixed values that are assigned to several regions globally [61]. In contrast, we used dynamic changes in the combustion ratio from daily FRP MODIS data. More fluctuations in C emissions from biomass burning were captured by FE in this study. As a result, the spatial distributions of FE in the study area were heterogeneous ( Figure 11) compared to the GPP and R e . Fire emissions are usually dependent on fire ignition events and are episodic [14]. For these reasons, relatively lower values were produced than the GPP and R e when described with the NBE (Figure 12). However, fires that destroy vegetation as potential carbon uptake sources and disturb the hydrological cycle through the subsidence of peatlands have the most impactful in changing the carbon balance of peatlands [62]. The FE in the Papua region was compared with GFED 4.1s ( Figure A1) because the visual interpretation of the map in the national scale was difficult. The FE of this study and GFED 4.1s in Papua were less detected than those of the other regions (e.g., Kalimantan). Relatively small scales fires were more detected in this study than by GFED 4.1s.
For future study, collecting more observations of carbon emission data during fire events is required for further validation. This study estimated FE based on the FRP of MODIS [46]; however, sub-daily scales of fire events might not be detectable. To estimate the total carbon emissions during one event from ignition to extinguishment, applications of geostationary satellite data are promising [46] to improve the FE due to their high temporal observation frequency.

NBE and Niño.3 Index
As mentioned above, R e and the GPP are generated by meteorological input variables, and datasets represented relatively homogenous distributions of 461.65 and 283.77 g C/m 2 /year for the spatial standard deviations when using every pixel, respectively. In the NBE equation, R e and the GPP compensated for each other. Therefore, both datasets hardly contributed to NBE fluctuations. The temporal changes in the NBE seemed to be dominated by variations in the fire emissions because NEE (R e − GPP) was overall constant (Figure 12). In addition, climatic impacts were analyzed by using the Niño.3 index, which represents how the sea surface temperature (SST) in 5 • N-5 • S 90-150 • E differs from a reference SST. The Niño.3 index provided by the Japan Meteorological Agency (JMA) showed El Niño or La Niña if the index value was positive or negative ( Figure 12). When the Niño.3 index was positive and negative, the averages of the NBE were 30.60 and 28.35 Mt C/month, respectively. FE in the El Niño and La Niña months were 2.58 and 1.25 Mt C/month, respectively. Usually, in the study area, during an El Niño period, the air temperature and precipitations were high and low, respectively-and vice versa in a La Niña period [63]. El Niño years are known for droughts [64] and more severe fire events than in normal years [65]. From the analysis of this study with the Niño.3 index, higher emissions were consistently found in the El Niño years. Furthermore, the other studies also revealed that higher CO 2 emissions occurred in El Niño years (2009-2010 and 2015-2016) from tropical regions through the analysis of the column-averaged CO 2 dry air mole fraction (XCO 2 ) by using the orbiting carbon observatory-2 (OCO-2), and the greenhouse gases observing satellite (GOSAT) [66][67][68]. Hence, higher R e values due to high temperature and more fire emissions were considered as the reasons for the high NBE in the El Niño periods.
However, the peak of the NBE tended to appear slightly ahead of Niño.3 peaks. Moreover, the peak of the NBE usually happened in the dry season from November to April, while the Niño.3 index was still increasing until April. Many fires occurred in the timing of on-set of the Niño.3 index. It might be concluded that the annual peak of the NBE appeared during the dry season of this region, and the timing of its peak was related to the on-set stage of the Niño.3 index.
In terms of R e , there were drastically decreasing periods: November 2009, March 2011, February 2016, and August 2019. The GWT data were investigated, since the GWT was only one factor for estimating R e in this study. The average (± standard deviation) of the GWT in the listed four months above and the other months were −0.14 (±0.08) and −0.25 (±0.34) cm, respectively. Anomalies of the two periods were also calculated. The averages of the anomalies of the selected periods (low R e ) (0.59 cm) were higher than the other months (0.26 cm). Therefore, we concluded that the low values R e in the listed four periods were induced by an abnormally shallow (high) GWT. Moreover, after the on-set of the Niño.3 index, the GWT was increased; that is, P was increased. Some previous studies have mentioned the impact of fire haze to generate more rainfall or fog in inducing cloud condensation nuclei [69,70]. Consistently, in this study, the low R e in the listed periods appeared in the month after huge fire emissions. In addition, Davison et al. (2004) [71] found that haze from fires could remain in the atmosphere for 20 days-to-one month after a fire event in Malaysia. Therefore, the low R e values happened in the month after large fire events, and the haze from the fires could have provided cloud condensation nuclei for one month after the fire events. Consequently, high precipitation associated with haze as resultant in a high GWT, which led low R e values.

NBE Reliability Comparing with Other Studies
The annual averages of the NBE were compared with estimations in other studies. Joosten [72] estimated C emissions from degraded Indonesian peatlands as 136.36 Mt C/year in 2008 by using an inventory-based method and reported a 345.27 Mt C/year NBE. In the same study by Joosten [72], a set of emission factors was applied to estimate global CO 2 emissions by multiplying the peatland areas for each country. Nation-wide comparisons were conducted in the study. However, the distribution of forest as the only source dataset hardly accurately depicts a carbon budget because SOM decomposition information could be omitted. Hooijer et al. [13] estimated 172.36 Mt C/year as a median of C emissions (with a possible range of 96.82-233. 19 Mt C/year) in Indonesian peatlands in 2006. The satellite-based net C budget (NBE) in this study estimated the value as 411.12 Mt C/year in 2006. In the same study by Hooijer et al. [13], the C emissions from peatlands were estimated by using a linear relationship between the GWT and CO 2 emissions from an empirically-derived look-up table that assigned CO 2 emissions to each land cover class. Though the estimation in this study was higher than in the other studies, this study has established a model to estimate various aspects including natural and anthropogenic emissions and absorptions, and this model was carefully validated with in situ observations by conducting spatiotemporal visualizations with satellite-based datasets.

Conclusions
In this study, we comprehensively evaluated the ecosystem carbon balance of tropical peatlands by using satellite-based data from 2002 to 2018 combined with in-situ observations in Indonesia. First, the GWT was estimated by using a typical drought index (KBDI). The implemented spatial distributions of the GWT were able to generate the spatiotemporal dynamics of R e . A MODIS-based FE dataset was used as a source of anthropogenic fire emissions. For the data for CO 2 uptake, a product of the MODIS gross primary production (GPP) was applied to calculate the NBE after modification with flux tower-based GPPs in Palangkaraya by using a robust regression method. We compared three types of flux tower observations for modeling R e , i.e., undrained forest (UF), drained forest (DF), and burned forest (BF). The highest accuracy among the three equations was found in the UF site. Both DF and BF showed a lower respiration than that of UF because the forests were already degraded and productivities had decreased. As such, we used a calibrated model that used UF to estimate the NBE time series. As a result of maps of various factors, R e and the GPP showed even distributions across the whole study area, although FE showed a heterogeneous distribution because the fires are episodic. The NBE was always greater than zero, meaning peatlands are a net emitter of carbon to the atmosphere. The average NBE was 353.51 Mt C/year from 2002 to 2018 when the peatland ecosystems were considered intact (without any disturbances). To investigate the impact of climatic changes, the Niño.3 index was compared with the NBE. El Niño periods, during which Niño.3 was greater than zero, might have accelerated R e and FE, which could have increased the NBE from tropical peatland ecosystems. The remote sensing-based diagnostic simulation of the C budget in this study provides spatiotemporal distributions of the net carbon balance, including natural emissions, anthropogenic-induced emissions, and natural absorptions in detail. The findings are expected to contribute to the creation of a feasible and robust assessment system for peatland ecosystems in terms of the daily carbon cycle by using satellite-based datasets. Table A1. Regression coefficient between a vegetation index (MODIS normalized vegetation index (NDVI)) and a biomass simulated by the vegetation integrative simulator for trace gases (VISIT) terrestrial biosphere model [28] in each land cover type of the MODIS land cover type product (MCD12Q1) to estimate biomass. Biomass = a * ln (NDVI) was used as a regression (where a is the regression coefficient). The fire emission estimation process was reported in detail in previous study [44].