Spatial Distribution of Forest Fire Emissions: A Case Study in Three Mexican Ecoregions

: This study shows a simpliﬁed approach for calculating emissions associated with forest ﬁres in Mexico, based on di ﬀ erent satellite observation products: the biomass, burnt area, emission factors, and burning e ﬃ ciency. Biomass loads were based on a Mexican biomass map, updated with the net primary productivity products. The burning e ﬃ ciency was estimated from a Random Forest Regression (RFR) model, which considered the fuel, weather and topographical conditions. The burned areas were the downloaded Maryland University MCD64c6 product. The emission factors were obtained from well-known estimations, corrected by a dedicated US Forest Service and Mexican campaign. The uncertainty was estimated from an integrative method. Our method was applied to a four-year period, 2011–2014, in three Mexican ecoregions. The total burned in the study region was 12,898 km 2 (about 4% of the area), producing 67.5 ( ± 20) Tg of CO 2 . Discrepancies of the land cover maps were found to be the main cause of a low correlation between our estimations and the Global Emission Database (GFED). The emissions were clearly associated to precipitation patterns. They mainly a ﬀ ected dry and tropical forests (almost 50% of all emissions). Six priority areas were identiﬁed, where prevention or mitigation measures must be implemented.


Introduction
The greenhouse gas (GHG) effect is a natural process that helps to maintain the thermal balance for life on Earth. However, the atmospheric composition is changing as a result of human activities, such as the burning of fossil fuels (coal, oil, natural gas), deforestation, and forest degradation. The consequent increase in GHGs has been associated with a rise in global temperatures [1][2][3], including in the period 2011-2015, which has been the warmest on record [4]. Wildfires contribute between 24% and 35% of carbon dioxide (CO 2 ), carbon monoxide (CO), Nitrogen Oxide (NOx) and methane (CH 4 ) emissions to the atmosphere, as well as a significant quantity of aerosols [5]. For this reason, and because of its role in vegetation changes, fire disturbance has been designated an Essential Climate Variable (ECV) [5] by the Global Climate Observing System (GCOS) to better understand Climate Change trends [6]. The GHGs associated with forest fires are CO 2 , CH 4 , and nitrogen oxides (NOx). CO is not a GHG, but it generates ozone (O 3 ) in the troposphere and has a deleterious effect on living organisms [7]. These gases are the result of chemical reactions resulting from the combination of three components: vegetation, heat, and oxygen [8]. Their emissions vary spatially and temporally, depending on the fuel moisture conditions: the more moisture, the lower the gas emissions (owing to incomplete combustion) [9,10].

Study Area
Three ecoregions in southwestern and southeastern Mexico, 22 • 0 -15 • 3 N and 86 • 30 -105 • 45 W, were selected for the study (Figure 1). The three ecoregions are forest dominated and susceptible to fires, but they differ in their natural characteristics and ecological response. These are 3 of the 22 ecoregions recognized by the ecoregion map of Mexico [23], and were selected according to the following criteria: (1) representative area of a forest ecoregion; (2) located within political states reporting a high incidence of forest fires (CONAFOR); (3) ecoregion with active fires identified with the Early Warning Forest Fires System, operated by the National Commission for Knowledge and Use of Biodiversity (CONABIO) [24]. The three ecoregions selected were: (a) Southern Sierra Madre (ecoregion 13.5), dominated by temperate forests (e.g., pine, oak, pine-oak, and cloud forests); (b) Southern Pacific Coastal Plain and Hills (ecoregion 14

Variables for Estimating Emissions
The selection of satellite and cartographic products for calculating the emissions was based on the factors proposed by Seiler and Crutzen (1980) and on Equation (1) [25]: , , = , , * , , * , * * 10 where , , is the amount of gas k released for a specific area (with the coordinates i, j) in Tg, , , is the biomass load for the same area (g m −2 ) (assuming the area has a homogeneous cover of fuel of vegetation type m), , , is the burning efficiency (0-1, dimensionless) of the fuel/vegetation type m, , is the burned surface of the same area (m 2 ), and (Emission Factor) is the amount of trace gas k released per unit of dry matter (g kg −1 of biomass). The data used for each parameter are described below.

Biomass
The goal was to establish a base map of the biomass and to update it for each year of interest. The distribution of the aboveground carbon stocks of woody vegetation in Mexico was selected as the base map [26,27], because it (1) uses field and remote sensing data with similar dates, (2) has been validated, and (3) covers the whole of Mexico with a 30 m spatial resolution. The biomass was estimated by applying conversion factors for the vegetation type, as proposed in [28], and it was resampled to 500 m.
Annual updating was performed using Net Primary Production (NPP), which is defined as the generation of organic material that plants accumulate in their tissue, as a result of photosynthesis by unit area [29,30]. NPP was extracted from MOD17A3 data [31] for the Gross/Net Primary Production (GPP)/(NPP). To estimate the biomass based on NPP, a conversion factor of 0.48 [28] was used. We

Variables for Estimating Emissions
The selection of satellite and cartographic products for calculating the emissions was based on the factors proposed by Seiler and Crutzen (1980) and on Equation (1) [25]: where M i,j,k is the amount of gas k released for a specific area (with the coordinates i, j) in Tg, BL i,j,m is the biomass load for the same area (g m −2 ) (assuming the area has a homogeneous cover of fuel of vegetation type m), BE i,j,m is the burning efficiency (0-1, dimensionless) of the fuel/vegetation type m, BS i,j is the burned surface of the same area (m 2 ), and E k (Emission Factor) is the amount of trace gas k released per unit of dry matter (g kg −1 of biomass). The data used for each parameter are described below.

Biomass
The goal was to establish a base map of the biomass and to update it for each year of interest. The distribution of the aboveground carbon stocks of woody vegetation in Mexico was selected as the base map [26,27], because it (1) uses field and remote sensing data with similar dates, (2) has been validated, and (3) covers the whole of Mexico with a 30 m spatial resolution. The biomass was estimated by applying conversion factors for the vegetation type, as proposed in [28], and it was resampled to 500 m.
Annual updating was performed using Net Primary Production (NPP), which is defined as the generation of organic material that plants accumulate in their tissue, as a result of photosynthesis by unit area [29,30]. NPP was extracted from MOD17A3 data [31] for the Gross/Net Primary Production (GPP)/(NPP). To estimate the biomass based on NPP, a conversion factor of 0.48 [28] was used. We applied the same conversion factor for the whole study area, because NPP is an annual product and has a coarse spatial resolution (1 km 2 ), while the biomass map was made for a specific period and has a higher spatial resolution. Both products, the base map and the annual biomass from NPP, were combined for an estimate of the biomass for a specific year.

Burning Efficiency
The burning efficiency corresponds to the fraction of the biomass or fuel consumed by the fire in a specific area; Mexican studies have used default values from general tables, or [19] have used information from the CONSUME model, which was designed specifically for the USA [19]. Here, we propose a Random Forest Regression (RFR) model, which has been extensively used in the last years for satellite classification, including forest fire studies [32]. This used the field data of fuel loads and their consumption, from 67 prescribed burnings in 2006 and 2007 (in 200 m 2 areas) in the following types of vegetation: oak, oak-pine, pine, pine-oak, cloud forest, shrubland, dry tropical forest, tropical forest, and tropical rainforest (data supplied by Dr. German Flores Garnica, National Institute of Forestry, Agriculture and Livestock Research [INIFAP]).
The burning efficiency was calculated with a Random Forest in the regression mode. The model was implemented in the R programming language. The variables used to create the model included fuel, weather conditions, and topography, all of which affect fire behavior. The following variables were considered: (a) the above-mentioned biomass carbon map, (b) a land use and vegetation map (1:250,000) produced by INEGI, (c) land cover maps for 2005 and 2010 produced with MODIS data at a 250 m spatial resolution [33], (d) Vegetation Continuous Field products (MOD44), downloaded from the MODIS web page, and (e) a time series of the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) at a 250 m spatial resolution in a composite of 16 days. We used the difference between the vegetation indices before and after each fire. In this case, both vegetation indices were used in order to select the best option. The weather conditions were obtained through a 100 h fuel moisture model, which used MODIS and Tropical Rainfall Measuring Mission (TRMM) data [34] for the temperature and precipitation and which applied the Tetens-Murray equation [35]. The model uses the flow of moisture into the types of fuels from one day to another [34]. The topography was expressed through a digital elevation model in terms of three indicators: altitude, slope, and aspect. The final RFR model consisted of 1000 trees. The individual tree predictions were then aggregated using the mean statistic. The coordinates of 67 areas of prescribed burning were characterized with 37 variables obtained from the above-mentioned cartographic and satellite products, in accordance with the date of the burning. The RFR model was assessed by means of the out-of-bag (OOB) estimations of accuracy (e.g., the correlation between the observed and predicted values), and the importance of the explanatory variables for the RFR model. In each run of the algorithm, the indicators were changed or removed depending on the importance of their participation in the RFR model. The method is described with more detail in [36]. The resulting RFR model was tested on the burned areas identified in March 2006.

Burned Areas
The burned areas were defined as the areas affected by forest or agriculture fires and are characterized by the presence of coal and ash, by the removal of vegetation, and by the alteration of the structure of the vegetation [6,37,38]. MCD45 and MCD64 satellite products, which are generated from the MODIS data, were evaluated and compared with the burned areas that were identified in high resolution images and polygons from fieldwork (CONAFOR). MCD64 was selected because it identifies burned areas better than MCD45 does; it displays burned and unburned pixels with a 500 m spatial resolution by showing persistent changes in the surface identified with the normalized difference between the MODIS bands 5 and 7, as well as by applying a dynamic threshold using the spatial and temporal distribution of active fires [39]. In October 2016, MCD64 Collection 6 was published with improvements. We validated this product with Landsat images from 2014 and 2017 by following the method described in [40].

Emission Factor
The emission factor is the amount of a compound released in grams per kilogram of fuel consumed (g/kg). The data in [41] have been widely used, but they are at the global level; thus, in this work, we used emission factors that were measured during a campaign in Mexico as part of the Megacity Initiative Local and Global Research Observation (MILAGRO) project, which was developed by the US Forest Service and Mexican agencies, including CONAFOR [42]. For this study, those results [42] were complemented with the work of [43], since the MILAGRO project did not consider two types of vegetation involved in the prescribed burnings.

Estimation of Emissions
Once the data were selected for estimating the parameters in Equation (1) The pixels identified as burned sites in MCD64A1 Collection 6 were characterized by the vegetation type, and the day of year (DOY) was transformed into the proper date. Pixels identified as natural vegetation were selected, eliminating those in cultivated areas. These selected pixels were characterized through a data cube generated with eight indicators of fuel and weather, according to the date of burning. The burning efficiency, which was calculated from the RFR model, was overlaid on the maps of the biomass and vegetation types; the latter is associated with the emission factors. The calculated emissions (Equation (1)) were integrated into cells of 10 × 10 km for their cartographic representation.

Uncertainty Calculation
The calculation of the four parameters to estimate the emissions with the remote sensing data entailed several assumptions that are considered to be valid; these assumptions were made because of a lack of field data, incomplete knowledge, and inherent data characteristics. Moreover, the spatial and temporal resolution and the minimum map unit generalize the information, so systematic and random errors are generated. Therefore, the uncertainty associated with the data can be estimated as [44]: where: U 2 i (i = 1, n) is the percentage uncertainty associated with each of the parameters, and U total is the percentage uncertainty of the product of the parameters.
The data used were (a) accuracy data of the biomass map [26]; (b) an accuracy matrix to validate the MCD64A1 product (burned area) using Landsat data, with a resulting accuracy of 0.79 and the calculation of the confidence interval of the producer accuracy [45]; (c) the prediction interval used to obtain the confidence interval, within an error of a critical value (t) and 95% confidence level; and (d) a confidence interval for the emission factors, calculated for each gas with an error formula with a critical value (t) and 95% confidence level.

Comparative Analysis
The results developed using the methods described above were compared with the GFED data. The GFED information has a spatial resolution of 0.25 • , and the percentage of burned area was estimated [46,47] to calculate emissions using the method proposed in [48]. The GFED annual files from 2011 to 2014, together with the GFED carbon emission algorithm, were downloaded. The algorithm is useful for calculating emissions by region and vegetation cover. The algorithm code was modified to generate a raster file of emissions within a cell of 0.25 • , since the original code only provides the results as a table. The carbon emissions were calculated from the measured data (annual biomass, burned area [MCD64 C.6], burning efficiency calculated with the RFR model), but the emission factors were based on the GFED information. The pixel results were incorporated into a cell of 0.25 • for a spatial comparison of the results.  Figure 3 shows the four annual biomass results, and Table 1 lists the correlation between the biomass map and biomass from NPP. Both biomass products have a similar behavior in each year, with differences in the quantity of the biomass as a response to meteorological conditions. The NPP data allowed for the interpretation of annual changes and the annual adjustment of the biomass variable.

Variables to Estimate Emissions
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 18 emission factors were based on the GFED information. The pixel results were incorporated into a cell of 0.25° for a spatial comparison of the results.  Figure 3 shows the four annual biomass results, and Table 1 lists the correlation between the biomass map and biomass from NPP. Both biomass products have a similar behavior in each year, with differences in the quantity of the biomass as a response to meteorological conditions. The NPP data allowed for the interpretation of annual changes and the annual adjustment of the biomass variable.

Burning Efficiency
For the Random Forest Regression (RFR) model, the OOB estimate yielded a correlation of 0.73 between the observed and predicted values, an absolute error of 0.079, and an RMSE (Root Mean Square Error) of 0.009. The model explained 49% of the variance in the data. The variables with a higher importance in the RFR model were (Figure 4): the vegetation type, land cover, biomass, VCF (Vegetation Continuous Field) tree cover, and 100 h fuel moisture; of these variables, the fuel moisture was the most dynamic variable (daily) and had a coarser spatial resolution. The least important variables were the VCF short vegetation cover, VCF bare ground cover and EVI; this last variable is dynamic, because it changes throughout the year.

Burning Efficiency
For the Random Forest Regression (RFR) model, the OOB estimate yielded a correlation of 0.73 between the observed and predicted values, an absolute error of 0.079, and an RMSE (Root Mean Square Error) of 0.009. The model explained 49% of the variance in the data. The variables with a higher importance in the RFR model were (Figure 4): the vegetation type, land cover, biomass, VCF (Vegetation Continuous Field) tree cover, and 100 h fuel moisture; of these variables, the fuel moisture was the most dynamic variable (daily) and had a coarser spatial resolution. The least important variables were the VCF short vegetation cover, VCF bare ground cover and EVI; this last variable is dynamic, because it changes throughout the year. Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 18 The variables had three different spatial-temporal characteristics: (a) Permanent or referential, such as the field data on the vegetation type, complemented with the map of the land use and vegetation, and the biomass data from a map [26]. (b) The annual variation, such as the land cover (2005), tree cover percentage, no tree cover percentage, and without vegetation cover percentage coming from VCF. (c) The intra-annual variation, which is associated with the relative differences in EVI of 16-day composites; the difference was calculated from one composite before and one after the fire (defined by the burned area date), and the 100 h fuel moisture for a specific date. These characteristics may be an influence on the burning efficiency behavior.
We used prediction intervals calculated by Quantile Regression Forest [49], as an evaluation element. The prediction intervals were the widest for tropical forests (maximum ±0.2), and the narrowest for temperate forests (minimum ±0.04). Temperate forests had the highest values of burning efficiency and a higher correlation than tropical forests. Table 2 lists the burning efficiency estimated by RFR model applied to burned areas of March, 2006, in contrast to the data presented by CONAFOR [50]. These estimates range from 0.37 in dry tropical forests to 0.61 in oak forests. The higher values correspond to oak vegetation, followed by tropical vegetation. The lowest estimates correspond to pine vegetation. There were not large variations between the estimates of the same vegetation type; in all cases the differences between the maximum and minimum values were smaller than 0.19. The values estimated by the RFR model were consistently lower than those of CONAFOR. The estimates most similar to those from CONAFOR were registered in tropical forests, and those that were most different were in dry tropical forests.

Burned Area
The improvements in MCD64 Collection 6 identified 52% more burned areas than version 5.1 did. The validation using Landsat imagery indicated an underestimation in the identification of the burned areas from the MCD64 product, especially when the burned areas were both few and small in extent, because of the coarse spatial resolution of the MODIS data. For instance, for the wet year 2014, an omission error of 80% was observed. The accuracy can be improved in years with more forest fires and large burned areas, such as in 2017, which has a 70% accuracy, with less than 30% in omission and commission errors. The variables had three different spatial-temporal characteristics: (a) Permanent or referential, such as the field data on the vegetation type, complemented with the map of the land use and vegetation, and the biomass data from a map [26]. (b) The annual variation, such as the land cover (2005), tree cover percentage, no tree cover percentage, and without vegetation cover percentage coming from VCF. (c) The intra-annual variation, which is associated with the relative differences in EVI of 16-day composites; the difference was calculated from one composite before and one after the fire (defined by the burned area date), and the 100 h fuel moisture for a specific date. These characteristics may be an influence on the burning efficiency behavior.
We used prediction intervals calculated by Quantile Regression Forest [49], as an evaluation element. The prediction intervals were the widest for tropical forests (maximum ±0.2), and the narrowest for temperate forests (minimum ±0.04). Temperate forests had the highest values of burning efficiency and a higher correlation than tropical forests. Table 2 lists the burning efficiency estimated by RFR model applied to burned areas of March, 2006, in contrast to the data presented by CONAFOR [50]. These estimates range from 0.37 in dry tropical forests to 0.61 in oak forests. The higher values correspond to oak vegetation, followed by tropical vegetation. The lowest estimates correspond to pine vegetation. There were not large variations between the estimates of the same vegetation type; in all cases the differences between the maximum and minimum values were smaller than 0.19. The values estimated by the RFR model were consistently lower than those of CONAFOR. The estimates most similar to those from CONAFOR were registered in tropical forests, and those that were most different were in dry tropical forests.

Burned Area
The improvements in MCD64 Collection 6 identified 52% more burned areas than version 5.1 did. The validation using Landsat imagery indicated an underestimation in the identification of the burned areas from the MCD64 product, especially when the burned areas were both few and small in extent, because of the coarse spatial resolution of the MODIS data. For instance, for the wet year 2014, an omission error of 80% was observed. The accuracy can be improved in years with more forest fires and large burned areas, such as in 2017, which has a 70% accuracy, with less than 30% in omission and commission errors. Figure 5 shows all of the burned areas identified by the MCD64 product in the period of 2011-2014. Spatially, zones with a concentration of burned areas were identified, as in the Guerrero and Oaxaca States into Southern Sierra Madre and Southern Pacific Coastal Plain and Hills ecoregions, or in the Yucatan and Campeche States into Plain and Hills of the Yucatan Peninsula ecoregion. Temporarily, there were large differences in the period; 2014 reported less than 1000 km 2 , in this case we expect fewer emissions. Meanwhile, other years reported more than 2000 km 2 of burned areas. The extreme case was in 2013, when were reported more than 5000 km 2 of burned areas, which could be translated into eight times more emissions that in 2014. Temporarily, there were large differences in the period; 2014 reported less than 1000 km 2 , in this case we expect fewer emissions. Meanwhile, other years reported more than 2000 km 2 of burned areas. The extreme case was in 2013, when were reported more than 5000 km 2 of burned areas, which could be translated into eight times more emissions that in 2014.

Emission Factors
The emission factors used in this work were selected for CO2, CO, CH4, NOx, and NH3 gases and PM2.5 particles. Table 3 lists these used emission factor. CO2 was the most abundant gas with values between 1603 g/kg in pine-oak forests and 1657 g/kg in tropical rainforests, and NH3 was the least abundant, with 0.5 g/kg in oak forests and 2.4 g/kg in tropical forests. The information was complemented with published data [43] for pine forests and desert shrubland.

Emission Factors
The emission factors used in this work were selected for CO 2 , CO, CH 4 , NOx, and NH 3 gases and PM 2.5 particles. Table 3 lists these used emission factor. CO 2 was the most abundant gas with values between 1603 g/kg in pine-oak forests and 1657 g/kg in tropical rainforests, and NH 3 was the least abundant, with 0.5 g/kg in oak forests and 2.4 g/kg in tropical forests. The information was complemented with published data [43] for pine forests and desert shrubland.

Spatial and Temporal Distribution of Emissions
The accumulated emissions of CO 2 , CO, CH 4 , NOx, and NH 3 gases and PM 2.5 particles for the four years in our study area were 67.5 (±20) Tg, over a surface of 12,898 km 2 . The uncertainty of emissions (29%) was lower than that reported elsewhere [19] (47%). Figure 6 shows the annual spatial distribution of emissions in cells of 100 km 2 , and Table 4

Spatial and Temporal Distribution of Emissions
The accumulated emissions of CO2, CO, CH4, NOx, and NH3 gases and PM2.5 particles for the four years in our study area were 67.5 (±20) Tg, over a surface of 12,898 km 2 . The uncertainty of emissions (29%) was lower than that reported elsewhere [19] (47%). Figure 6 shows the annual spatial distribution of emissions in cells of 100 km 2 , and Table 4 lists the total emissions.    The emissions were highest in 2013 (27.5 Tg), accounting for 40.8% of the four-year total, with a burned area of 5272.75 km 2 . The spatial distribution was more uniform in 2013 than in the other years ( Figure 6C): 53.8% in the Southern Sierra Madre ecoregion (13.5); 28.1% in Plain and Hills of the Yucatan Peninsula (15.2); and 18.1% in Southern Pacific Coastal Plain and Hills (14.5). In contrast, in 2014, the emissions were 2.9 Tg, representing 4.2% of the four-year total, with a burned area of 657.25 km 2 . Furthermore, the concentrations of gases were low ( Figure 6D)

Emissions and Precipitation
Emissions were associated with the distribution of precipitation in Mexico [53] (Figure 7). Two seasons were identified: wet from May to October; and dry from November to April.

Emissions and Precipitation
Emissions were associated with the distribution of precipitation in Mexico [53] (Figure 7). Two seasons were identified: wet from May to October; and dry from November to April. The concentration of emissions in the dry season has two main causes: (a) the vegetation's susceptibility to burning in the driest months, and (b) the use of fire in agricultural activities. In Mexico, the months with the highest activity of agricultural burning are those in the early part of the year (e.g., in 2017 the prescribed burning in Campeche State was from 12 March to 31 May).
The month with the highest emissions was one month after the driest month and before the wet season; therefore, May, April, and March had >90% of all emissions. The remaining 10% was distributed mainly among the other months of the dry season, with a smaller proportion in the wet season months. In none of the four years were there emissions recorded in September, the month before the end of the wet season and often the month with the highest precipitation. Only two years showed any emissions in August (0.03%) or October (0.01%). The 2013 wet season was the wettest of the four years; hence, in 2014, fires were fewer and the emissions were the lowest (Figure 7). The dry season in 2013 was the driest, and it had the highest emissions.    The month with the highest emissions was one month after the driest month and before the wet season; therefore, May, April, and March had >90% of all emissions. The remaining 10% was distributed mainly among the other months of the dry season, with a smaller proportion in the wet season months. In none of the four years were there emissions recorded in September, the month before the end of the wet season and often the month with the highest precipitation. Only two years showed any emissions in August (0.03%) or October (0.01%). The 2013 wet season was the wettest of the four years; hence, in 2014, fires were fewer and the emissions were the lowest (Figure 7). The dry season in 2013 was the driest, and it had the highest emissions.   (Figure 1). Dry tropical forests have lower emissions but a greater variability, from 8.2% to 17.4%. In 2012, emissions were distributed more uniformly in temperate forests, such as pine-oak, oak-pine, and oak in the 13.5 Southern Sierra Madre ecoregion, with the percentage in tropical forests at <15%. In cloud forests, which are restricted to the mountainous area in that ecoregion, the burned pixels ranged between 2.7% and 4.9% over the four years. Meanwhile, the emissions generated in tropical rainforests and hydrophilic vegetation represented less than 1% for all of the years; the former was restricted to small areas of the ecoregions 15.2 and 13.2, and the latter was limited to the coastal zone of the ecoregions 15.2 and 14.5. Grassland produced emissions at low percentages, between 3.7% and 5.9%, in all periods with the exception of 2014, when they reached 12.0%.

Emissions and Vegetation Type
The annual emission dynamic depends on the distribution of vegetation and precipitation. In January, the vegetation type that had more emissions was grassland, because it is light fuel and quickly loses the moisture that it accumulates during the wet season. As each year progressed, emissions significantly increased in temperate forests, reaching more than 50% in May and June, while in tropical forests the emissions decreased. In July, emissions were generated in the tropical forests. By the end of the year, emissions had again increased in temperate forests and grassland. Therefore, under the conditions described here, the distribution of emissions at the ecoregion level is due mainly to environmental characteristics.

Areas of Interest in Terms of Emissions
Six areas were considered to be of interest because of a concentration of cells with emissions exceeding 119 Gg and because they were spatially continuous ( Table 6 and Figure 8). Three areas were in temperate forests and the remainder were in tropical forests. The origin of forest fires in these areas is associated with agricultural activities and urban expansion; e.g., in Cancun, the cells with a concentration of emissions were along the highway and close to the urban area. When fire used in agricultural and urban activities loses control and involves natural vegetation, it is considered a forest fire.

Comparative Analysis of Results with Global Fire Emissions Database (GFED).
The present estimates of carbon emissions during 2011 and 2013 were similar to the GFED data for forested areas (Table 7). However, in all years except 2011, there was a lack of correlation because of strong local discrepancies. The highest coincidence in the estimated amount of carbon was in 2013 in Plain and Hills of the Yucatan Peninsula (ecoregion 15.2), and the lowest was in 2012 in Southern Sierra Madre (13.5). In the latter, the cells with the highest emissions identified by GFED were towards the south of the cells with the highest emissions identified by this work. The MCD12Q1 land cover product, which is an input of GFED, identifies this zone as woody savannas, whereas the INEGI cartography of Mexico identifies it as oak forests, located on the leeward side of Southern Sierra Madre, where the climate is drier than that of the windward side. In the GFED, the savanna is not considered to be a forest. Hence, differences in the cartography of the Southern Sierra Madre are the main source of discrepancy in the results. In contrast, the correlation was higher when there were more emissions in the Plain and Hills of the Yucatan Peninsula (15.2) ecoregion. Another element to consider is the difference between the versions of the burned area product; the present work used the latest version (Collection 6), whereas GFED uses a previous version.

Comparative Analysis of Results with Global Fire Emissions Database (GFED)
The present estimates of carbon emissions during 2011 and 2013 were similar to the GFED data for forested areas (Table 7). However, in all years except 2011, there was a lack of correlation because of strong local discrepancies. The highest coincidence in the estimated amount of carbon was in 2013 in Plain and Hills of the Yucatan Peninsula (ecoregion 15.2), and the lowest was in 2012 in Southern Sierra Madre (13.5). In the latter, the cells with the highest emissions identified by GFED were towards the south of the cells with the highest emissions identified by this work. The MCD12Q1 land cover product, which is an input of GFED, identifies this zone as woody savannas, whereas the INEGI cartography of Mexico identifies it as oak forests, located on the leeward side of Southern Sierra Madre, where the climate is drier than that of the windward side. In the GFED, the savanna is not considered to be a forest. Hence, differences in the cartography of the Southern Sierra Madre are the main source of discrepancy in the results. In contrast, the correlation was higher when there were more emissions in the Plain and Hills of the Yucatan Peninsula (15.2) ecoregion. Another element to consider is the difference between the versions of the burned area product; the present work used the latest version (Collection 6), whereas GFED uses a previous version.

Discussion
The use of satellite Earth observation data allowed for the spatial and temporal estimation of emissions, providing an annual updating of the biomass and burned areas, as well as the evaluation of environmental conditions for estimating the burning efficiency.
The parameter with the highest uncertainty was the biomass load (27%). Even though the biomass map had an RMSE of 54% in its validation [26], it was still considered acceptable because similar results have been reported by other projects; e.g., the GlobBiomasss project, University of Leicester (2017), reported RMSE values of 51%, 53% and 55% for the years 2005, 2010, and 2015, respectively, in the Central Mexico and Yucatan Peninsula regions. As a complement to the present work, the emissions were estimated by using only the biomass map [26] for each of the four years; this led to a value that was lower by 19 Tg, thereby reinforcing the importance of annual updating.
The burning efficiency results obtained for March 2006 were within those values reported elsewhere: (a) in the literature [50]; (b) in the documentation of prescribed fires [54], and (c) in a study to estimate the spatial variability in burning efficiency in Spain [55]. Hence, they were considered appropriate, although they were lower than those reported in the literature. The results also indicated that the behavior of the burning efficiency can be explained by the vegetation type; temperate forests have a relatively low species diversity, and the fire behavior is more uniform than in the species-diverse tropical forests. Here, field data are a fundamental aspect of the method. We considered that the Random Forest approach [36] offers an approximation for estimating the burning efficiency over a large and biodiverse territory such as Mexico; it considers spatial and temporal changes in environmental conditions when a fire occurs, without depending on a unique and static value.
Despite the relatively low accuracy revealed by the omission error for the burned area product, MCD64 has been the best operational product with the historical and recent data, because the data are processed with the same algorithm. This is supported by a study [56] that concluded that the Collection 5.2 MCD64 was the best of three evaluated burned area products, and Collection 6 identifies more burned areas than Collection 5.2.
The emission factors obtained from the MILAGRO project are higher than those reported in the global study of [41] for CO 2 and CH 4 but lower for CO. These differences have implications in the estimation of emissions because CO 2 is the gas with the highest emission factor. However, the MILAGRO data are specific to Mexican vegetation.
Southern Sierra Madre (ecoregion 13.5) had the highest percentage of emissions, but the vegetation type most affected was the tropical forest distributed in ecoregions 15.2 (Plain and Hills of the Yucatan Peninsula) and 14.5 (Southern Pacific Coastal Plain and Hills). If the emissions originating in temperate forests are added, the result is similar in tropical forests. Therefore, it is important to highlight two elements in the calculation of emissions: the biomass quantity and burning efficiency. Although the amount of emissions per ecoregion reflects the area of that region, it is also influenced by the type of vegetation, the date, and the causes of the fire. For example, ecoregion 13.5 represented 34.8% of the study area and 54.8% of the emissions over the period, whereas ecoregion 15.2 represented 44.3% of the study area but only 28.7% of the emissions.
The differences found between the GFED data and our results are due to the scale of the analysis; GFED is at a global scale without regional detail, giving a general view of the problem, whereas our results have a local precision. We also identified an important difference between our results and a number of previous works [20][21][22]. The results of those three works are similar to each other (9.6 Tg, on average) but lower than those of this study. We found that the burned area parameter had a large influence. CONAFOR statistics recorded less than twice the burned area identified by MCD64A1 Collection 6.0, and this can be explained by the discrepancy in the definition of forest fires. Many fires involving vegetation are associated with the clearing of agricultural fields and with forestry operations; these are not considered uncontrolled fires, and hence are not registered as forest fires, but they are still recorded by MCD64. Hence, it is necessary to identify whether the emissions come from a forest fire or from a change in the land use [57]. Another parameter with an influence on the differences was the biomass. When only the biomass map was considered, without the annual update of NPP, the calculated mean of the emissions fell from 16.8 Tg to 12.1 Tg; even so, this is still a high value.
The spatial and temporal distributions of emissions in the study area were associated with the distribution of precipitation, drought, and human activities. There were two very important elements in the emissions distribution: hydrological drought and biomass. The values of the estimated emissions were expected to be low because the estimates of the burning efficiency were small compared with those of other studies. However, the emissions were higher, owing to the inclusion of different inputs, such as a yearly update of the biomass and burned areas. Therefore, it is important to consider the changes over time in the variables and the contribution of satellite data, e.g., the observation of the phenomenon on the day of occurrence. In this work, we considered drought conditions by using the fuel moisture model and the annual update of the biomass.

Conclusions
The availability of satellite data and products allows for the estimation and update, both dynamically and periodically, of three out of four parameters for calculating the emissions from forest fires. However, there is still a lack of data with the same resolution and scale to contribute to the precision of the spatial analysis. This was illustrated by a comparison between the GFED data and our results, given that the values of the emissions were very similar, although local discrepancies were observed when data with a better spatial resolution were used.
The present study contributes a product of burned areas that identifies areas affected by fires in vegetation, irrespective of the forest fires registered in the statistics. It is important to collaborate with CONAFOR in using remote sensing data to improve the statistics. Recently, CONAFOR has improved the method for estimating burned areas in the field, and has sometimes incorporated satellite data, in collaboration with CONABIO. Remote sensing data and field knowledge will increase the quality of information regarding emissions. We demonstrated the use of remote sensing data and machine learning to calculate parameters such as the burning efficiency, taking into account natural conditions on the date of the fire; this allows for a better understanding of emissions resulting from forest fires. The areas of interest with burned biomass emissions were located at a larger scale than that used in other works and included an adjustment to the environmental conditions on the date of the fire. This information can help to implement measures to prevent or to mitigate fires according to the environmental characteristics.
The described methodology may be implemented on a countrywide scale because the variables and field data that it uses in the Random Forest Regression (RFR) model are to some degree flexible. As mentioned, it includes an evaluation of the accuracy of its estimates. This allows the methodology to be applied on other ecoregions. Finally, the satellite imagery proposed as input is generated periodically and is available on different distribution platforms.