Assessment of carbon dynamics in Ecuadorian forests through the Mathematical Spatial Model of Global Carbon Cycle and the Normalized Differential Vegetation Index (NDVI)

Conservation and sustainable development of forests are mitigation mechanisms against climate change due to the forest carbon sink capacity. Therefore, biomass estimation allows to assess forest productivity and control carbon budgets. In Ecuador, biomass and carbon sequestration studies are scarce. Thus, we estimated and forecasted changes in biomass of Ecuadorian forests through the Mathematical Spatial Model of Global Carbon Cycle and the Normalized Differential Vegetation Index. The mathematical model describes the processes of growth and decay of vegetation in terms of carbon exchange between the atmosphere, plants and soil under anthropogenic impacts. The vegetation map and the biomass of 2017 (4,86 Gt) were developed with remote sensing methodology in ENVI 5.3 and ArcGIS 10.3 programs. The observed biomass decrease between 2000 and 2010 was due to the high deforestation rate. Thanks to conservation and reforestation policies and the compensatory effect between the atmosphere and forests, a biomass increase is expected until 2060. According to the vegetation map, Amazon region has a better plant vigor, followed by Andean and Coast regions, where scattered vegetation predominates. This information is useful for planning environmental practices such as forest conservation and reforestation in order to increase carbon storage.


Introduction
Global climate change is the result of the imbalance between population increase and the natural resources ability to sustain the growing demand. As a result, anthropogenic impacts such as deforestation, soil degradation and atmosphere pollution have increased the levels of CO2 and triggered the global warming. In Ecuador, CO2 emissions from transport and industries as well as deforestation and soil erosion have contributed to the increase of this gas. According to the Food and Agriculture Organization of the United Nations, Ecuador had one of the highest annual rates of deforestation (1.8%) between 2001 and 2010.
Conservation and sustainable development of forests are mitigation mechanisms against the adverse effects of climate change. Forests act as carbon sinks, absorbing and storing C in biomass (organic matter) and soils. Approximately, the 50% of CO2, captured by photosynthesis, is fixed as organic carbon in trunks, branches, leaves and roots. Therefore, the biomass estimation allows to assess the forest productivity, control carbon budgets and have a better understanding of carbon cycle [1].
The assessment and forecast of carbon dynamics and its biospheric stability are scarce and non-existent in Ecuador. Just few biomass studies have been developed about forest and the biomass in 2012 is detailed in the National Forest report submitted to the FAO in 2015 [2].Therefore, the objective of this research is to estimate and forecast changes in the ecosystem function of carbon sequestration in Ecuadorian forests and their balance under anthropogenic impacts through the biomass estimation.
Remote sensing and global carbon models are indirect methods to estimate biomass. Carbon cycle model is the only method that forecast future variations. In comparison with the direct methods they are not invasive and costly. Besides, for estimate biomass in large forest areas they are not time and labour consuming. Therefore, remote sensing and the model of global carbon cycle are used in this research as a better methodology to forecast the carbon dynamics and know the role of terrestrial ecosystems in anthropogenic compensation [3].

Remote sensing and Normalized Difference Vegetation Index (NDVI)
Earth Remote Sensing (ERS) is a non-physical contact method of acquiring, processing and interpreting images of an object taking from distance by aircrafts and satellites. Some of the uses of ERS are the assessment of the state of the territory with Vegetation Index (VI) and the quantification of some properties of the earth's surface like the biomass. Studies have applied VIs to estimate biomass. For instance, Dong et al. [4] correlated the vegetation greenness of temperate zones with biomass field inventories for long time series. One of the VI widely used is the NDVI because it minimizes the topographic effects while produces a linear measurement scale. In addition, division by zero errors are significantly reduced. NDVI was introduced by Rouse et al. [5] and takes advantage of the fact that greener or healthier vegetation absorbs more visible light and reflects a large amount of near infrared light, while unhealthy or sparse (less green) vegetation reflects a large portion of visible light and less near infrared light. The NDVI formula is [6]: NDVI = (NIR-RED)/(NIR+RED) Where NIR is the reflectivity in the near infrared zone (0,75 -1,3 μm) and RED is the red zone (0,62 -0,75 μ m). The relations between the coefficients in the formula generates values in the range of -1 and +1, where the negative values represent bodies of water, values slightly greater than 0 belong to bare ground, values between 0,2 to 0,8 represent spread vegetation, greater than 0,8 correspond to dense green vegetation and those values close to 1 are for rainforests. The index depends on the vegetation composition, its proximity, exposure and the inclination angle of the surface and can vary depending on the use of the soil, phenological season, water situation of the territory and the weather [7].
One limitation of NDVI is the saturation in cases of dense and multi-layered canopy. Also, it cannot minimize the effects of soil bottom; it means that a certain proportion of their values, negative or positive, represents the brightness of the soil. The lighting is another drawback because the atmospheric changes do not affect all wavelengths in the same way [8].

Spatial Model of the Global Carbon Cycle Model
Carbon cycle is the carbon fluxes exchanged between four deposits: the atmosphere, the biosphere, the oceans and sediments. This cycle has maintained the carbon balance on the Earth and the stable temperature. However, human beings have disturbed the equilibrium of this. This unbalance is the main reason for the increasing interest on the Carbon Cycle [9].
The modelling of Carbon cycle allows to forecast the future behaviour of carbon in the climate system and look for sustainable development. Models like the spatial model of global carbon cycle in the Atmosphere -Plants -Soil (APS) system are focused either on detailed description of elements of biogeochemical cycles or on investigation of new data of measurements.

Normalized Differential Vegetation Index
In total 17 Landsat 8 images were downloaded from Earth Explorer of the United States Geological Survey (USGS) in 2017. The clearest images and with least visual interference over the vegetation were visually selected. The processing of Landsat images (radiometric calibration and atmospheric correction) in the ENVI 5.3 program required the multispectral file (MTL).
The NDVI was calculated in ENVI 5.3. A scale of 5 intervals was assigned to classify the vegetation: less than 0 includes no vegetation areas like bodies of water, flooded areas and clouds, from 0 to 0,1 corresponds to bare soil, degraded land, roads, settlements, soil without vegetation cover, from 0,1 to 0,3 is scattered vegetation, grass, scattered shrub, irrigated crops and plow fields), from 0,3 to 0,5 is open vegetation (forest plantations, bushes, slow-growing plantations and dry forest) and greater than 0,5 is closed vegetation (dense growth plants, evergreen forest and montane forest). NDVI ranks were based on those proposed by Cartaya et al. [10]. The final NDVI map of Ecuador was created in ARGIS.

Quantitative calculation of phytomass
The forest biomass estimation is based in the methodology applied by Grigorets et al. [3]. The reference biomass per unit area for 2012 was 0,0283 t/m2 [2]. With this value the amount of phytomass for a pixel of 900 m2 (resolution of Landsat satellite is 30 m) was 25,5 t. The value of 25,5 t corresponds to 1 on the NDVI scale, which means the optimum value of vegetation status. The phytomass corresponding to each pixel and according to its value of the NDVI was calculated by multiplying the number of pixels by its NDVI, established in the NDVI histogram ( Fig. 1) and by 25,5 t. The process was repeated for all the Landsat images. Finally, the phytomass values of the 17 scenes were added, giving a total of 4863900904 t for 2017 year.  , C is total carbon, phytomass is and the soil humus is . Carbon dioxide emissions are represented by , the amount of carbon in the decomposed biomass is and is the amount of carbon released by soil erosion.

Assessment and forecasting of phytomass and humus dynamics
You The assessment and forecasting was developed under the Spatial Model of the Global Carbon Cycle. It is a universal model (for all ecosystem) and divides the Earth into cells of size 0.5 x 0.5 degrees in the geographic grid (50 x 50 km). The model describes the processes of vegetation growth, decay, accumulation and humus decomposition in terms of carbon exchange between the atmosphere, plants and soil in every cell (Fig. 2). It contains more than 100 thousand nonlinear differential equations.
In each land cell, the ecosystem is characterized by the carbon amount per area that exists in the phytomass ( ) and in the soil humus ( ). The area of i is denoted by σi and C represents the total carbon in the atmosphere. The climate is characterized by the annual surface temperature and the annual precipitation , they are calculated using a spatial climate model of the general circulation in the atmosphere and ocean. It is assumed that the annual vegetation production in each cell depends on C, as well as and and does not depend on the type of ecosystem. The time unit is one year. The rate of decomposition of humus, is represented as a function of and , . The model takes into account the three anthropogenic factors that increase the CO 2 concentration. The first is the carbon dioxide emissions generated by fossil fuels burning . The second factor is deforestation. The amount of carbon in the decomposed biomass is denoted by . The third factor is soil erosion. The dead organic matter is decomposed with a certain delay and goes to the atmosphere through the ocean; is the amount of carbon released by soil erosion. Deforestation and erosion data were taken from the Environment Ministry of Ecuador from 2000 to 2013.
The carbon dynamics in the Atmosphere-plant-Soil model is explained in the following equations: (1) Where is the annual production of plants per unit area in the cell , and are coefficients. There are two expressions that depend on the primary annual production network. In the first, annual production depends on the concentration of CO 2 in the atmosphere and the phytomass of the vegetation.
(2) Where is the annual production of the initial state, is the phytomass of the vegetation expressed in carbon units. and are initial values corresponding to the respective variables and β is a coefficient. The second expression was developed by Tarko, where the annual production depends on the atmospheric concentration of CO 2 and not on the phytomass.
(3) The function expresses non-linear dependence on annual production as a function of temperature and precipitation. The model has been complemented by the carbon cycle model in the atmosphere -ocean system. It is considered in the model that the annual CO 2 emissions of fossil fuels are mixed in a latitudinal direction for approximately two weeks and from 2 to 3 months in a southerly direction. Therefore, the territory is subject to the effects of climate change, depending on the total emissions for one year. To calculate the climate change generated by global warming, the country's total emissions and the special climate model of the entire planet were used [11].

Results and Discussion
The NDVI map (Fig. 3) shows the current state of vegetation. In general, coast region is classified as bare land, degraded land and scattered vegetation. The agriculture expansion in the coast, especially with banana and cacao crops, between 1990 and 2000 is the main reason to the current state. In addition, the shrimp industry has eliminated the 70% of mangroves. Also, African palm crops have increased in the last decades in the north coast. The same scenario is observed in the north Amazon due to oil and wood extraction. Nevertheless, the majority of Amazon region is classified as closed vegetation. As a result an important carbon sink is located here. The Amazon forests conservation is because approximately 40% of the amazon is National Protected Area. Finally, in the Andean region there is a mixture of vegetation states. Although, agriculture and livestock predominates in this area, there are areas with a NDVI more than 0.5 owing to its difficult access (altitudes more than 3000 m.a.s.l).  The correlation between the results (Fig. 3) and the spatial distribution of average biomass values developed by Anaya et al. [12] is high. They established the Colombian aerial biomass with MODIS images and in the final map the coast got the least values of biomass (0 and 55 t/ha), followed by Andean region and finally the Amazon with the highest values of 500 t/ha. Moreover, the Ecuadorian NDVI classification matches with Meneses [13] analysis where the greatest NDVI values correspond to high and medium altitude tropical forests and the lowest NDVI values for deserts, pastures and bushes. Also, Meneses found that the NDVI values change according to the season (with precipitation or dry period).
Carbon storage strongly depends of the forest density. Dense forests have a greater capacity to store carbon than open forests [9]. So, a better efficiency for accumulating carbon in Amazon and Andean regions is expected than in the coastal forests. The factors that affect the absorption and storage of carbon in the forests are the temperature, precipitation, mass density, soil, slope, height, topographic conditions, growth index and age [9]. Those environmental factors differ in the three natural regions of Ecuador: Coast, Sierra and Amazon regions (Fig. 3). In the coast region the altitude varies from 0-800 to 1000 m.a.s.l (meters above sea level) and the climate is warm, cool and dry in the center and south and hot and humid in the north. The Sierra or Andean region has different climatic zones: the Andean tropical with average temperature between 20 and 25 °C and 1500 m of altitude. The Andean subtropical from 1500 to 2500 m.a.s.l with average temperature of 20 °C. The temperate zone that reaches up to 3500 m.a.s.l with temperature of 17 °C. The cold floor from 3 500 to 5 650 m whit a variable temperature between 1 and 10 °C and finally the glacial floor with temperatures below 0 °C. The Amazon has a hot humid climate with temperatures between 22 and 26 °C and abundant rainfall (more than 3000 mm per year) and its altitude reaches 800 m.a.s.l. Campo et al. [14] affirm that high rate of average annual precipitation (AAP) favors the growing of tree, hence, the highest biomass and carbon sequestration rate in Amazon forest are the result of the fast grown due to warm and wet climates with long growing seasons [9]. On the other hand, the lower carbon storage capacity in coastal forests is attributed to low precipitation.
The lower amount of humus and phytomass estimated by the model between 2000 and 2010 (Fig. 4) is attributed to the past high rate of deforestation of 1.8%. A slight increase appears between 2010 and 2020, a possible explanation is the deforestation decrease in recent years. In fact, between 2014 and 2016 the deforestation rate (61 thousand ha / year) has diminished in 33,7 % in comparison to the deforestation rate of 1990 -2000 period (92 thousand ha / year). This advance in the forest conservation is due to the current government policies that consider the nature as a subject of law in the modified constitution of 2008. Moreover, Ecuador is a member country of the REDD program (Reduction of Emissions from Deforestation and Forest Degradation in Developing Countries) since 2011. In addition, Socio Bosque Program (PSB) conserves forests through agreements with private forests and indigenous communities. It would be expected that these programs will be the main reasons for the future phytomass and humus increase (6.06%) towards 2060.
The compensatory effect between the atmosphere and the forest ecosystems also contributes to the humus and phytomass increase. The effect suggests a greater production and growth of the plant phytomass when the CO 2 emissions increase. Studies with woody plants exposed to high concentrations of CO 2 corroborate this assumption of compensatory effect. The trees showed a greater and faster growth with a bigger leaf area, higher rate of photosynthesis and more organic matter in soil [15]. Therefore, it is expected that the biomass of Ecuadorian forests increase with the raise of carbon dioxide emissions. Similarly, the same trend has been reported for countries like Brazil, China, India, Russia and African countries during 2000-2060 [16].  The forecast agrees with the results of carbon fluxes dynamics of tropical lowland forests and tropical montane forests in South Ecuador. Those studies applied the forest growth model FORMID with the strength on spatial and temporal extrapolation. The forest acted as a carbon sink with a maximum net ecosystem exchange of 9,3 Mg C (ha/year) during its early successional stage between 0 to 100 years [17].
The biomass estimation with the Global Carbon cycle showed a high correlation with the reference biomass and the calculated with NDVI. The biomass in the model is 4,4674 Gt for 2010 and the reference biomass is 3,22 Gt for 2012. For the year 2020 it is 4.4686 Gt and the estimated with the NDVI is 4.86 Gt for 2017. It is clear in both methods the biomass increase over the time.

Conclusions
The lower amount of humus and phytomass in the carbon dynamic of Ecuadorian forests from 2000 to 2010 is consistent with the rate of deforestation and erosion. In the following 10 years a slight increase in biomass is observed as a result of two possible reasons 1) Conservation policies and 2) Compensatory effect between the atmosphere and forest ecosystems. Finally, the trend for 2060 is an exponential increase in biomass, exceeding the phytomass values of 2000 in 6.06%. This increase is also attributed to the same two reasons.
Both methodology, mathematical model and remote sensing, showed a high correlation in their values It increases the confidence of the model predictions. In addition, those methods do not affect the ecosystem and reduce costs and time. However, it is recommended to combine with field verification to optimize the results.
A carbon sink is concentrated in the Amazon region with NDVI greater than 0.5, followed by the Sierra. On the other hand, the Coast is classified as bare soil (NDVI 0 -0.1) and scattered vegetation (NDVI 0.1 -0.3) due to the past agriculture expansion.
A theory of forestry transformation is identified thanks to the implementation of conservation and reforestation projects since 2008, such as PSB and REDD+. It implies an increase in forest cover and, consequently, an increase in carbon sinks through the sustainable forest management. Also, the inclusion of the rights of nature in the Ecuadorian Constitution has allowed a better forest management.
The reforestation strategy must consider the type of vegetation and the climatic floor in which it is developed to avoid release CO 2 instead storage it. In the past to reduce the concentration of atmospheric CO 2 wrong solutions were applied like reforestation with exotic species such as eucalyptus and pine. This led to a hydrologic unbalance that increased the soil aridity and diminished the carbon sequestration potential.
The biomass estimation and its forecasting can be used in the planning of environmental practices such as REDD+ reforestation activities, national strategies for the forest conservation, as maps of recoverable forests with different possibilities of carbon capture in different geographical and climatic conditions, and also in the framework of the implementation of the concept of "avoided deforestation".