Interactive comment on “ Simulating Carbon Sequestration using Cellular Automata and land use assessment ; Case of : Karaj City , Iran ”

First of all, we would like to thanks for very great comments of reviewer (1) to help us improve the manuscript and to inform you that the revised version of our manuscript with considering all comments of reviewer (1) has been prepared and it is ready for submission to your journal. Based on comments, the outline of corrections is as follows: You should improve the abstract, see line 18 and 19 you can rephrase these lines. The abstract revised based on your comment. It starts out with one line introduction, the study objective, the methodology, main finding and conclusion. Line 18 and 19 revised. These changes have been highlighted in the text.


Introduction
Carbon sequestration is defined as the process of removing carbon from the atmosphere and depositing it in a carbon reservoir (UNFCC, 2015).Soil carbon sequestration is the process of storing carbon in soil, sometimes for thousands of years.When carbon is removed from the atmosphere and stored in the soil as organic carbon, its contribution to global warming is reduced.Plants absorb carbon from the atmosphere through photosynthesis and they release carbon to the atmosphere through respiration.Any carbon that is stored in plant tissue is either consumed by animals or, after the plant dies, it enters the soil as decomposed organic matter.Soil organic matter is a complex mixture of carbon-containing materials that include decomposed plant and animal tissues, microbes (protozoa, nematodes, fungi and bacteria) and carboncontaining minerals.Factors such as climate, natural vegetation, soil texture and drainage all influence the level and duration of soil carbon sequestration (Schlesinger, 1984).
Humic organic materials form the main sequestered constituents in soil, whereafter humus is transformed into organic and inorganic compounds that become located deep down in the soil due to agricultural plowing, root transportation by plants and calcareous soil.Carbon sequestration has been proposed as a way to slow the atmospheric and marine accumulation of greenhouse gases that are released by burn-ing fossil fuels (Hodrien, 2008).Terrestrial carbon sequestration is the result of a balance between the different stages of the carbon cycle in the biosphere and pedosphere, such as photosynthesis, plant growth, congestion and carbon accumulation in soil, as well as carbon emissions from breathing organisms, microbial decomposition of leaf litter, and oxidation of organic carbon in soil and land degradation.Several factors are involved in this process, which can be classified into two main classes: physical and managerial.Physical factors can be further divided into three subclasses: soil, climate and terrain geometry.Since the physical factors are generally difficult, sometimes impossible, to control, management of carbon sequestration must focus on managerial factors (Ferreiro et al., 2010).
Because the majority of carbon is sequestered in soil, organic carbon that is contained in soil plays an important role in the overall sequestration process.The amount of carbon that is sequestered in soil changes significantly according to location, topography, bedrock or vegetation type and previous management approaches (Schoonover and Crim, 2015).Furthermore, the amount of carbon found in the roots, litter and biomass of microorganisms in the soil varies according to the time of year due to growing season effects and time required for biomass decomposition.Erosion, compaction and reduction of soil permeability, loss of soil structure, mineralization and oxidation of humic substances all contribute to reducing the amount of organic matter in soil (Bruce et al., 1999).Although the amount and rate of carbon sequestration is higher in mild and humid tropical forest ecosystems, the high level of plant respiration that is associated with high humidity and temperature results in a fast release of carbon dioxide back into the atmosphere so that, on balance, the net carbon storage is reduced.Therefore, arid and semi-arid regions are considered to be the most desirable locations for carbon sequestration (Lal, 2008).Indeed, international organizations, such as FAO and UNDP, have chosen to implement programs of carbon sequestration in arid and semi-arid regions in an effort to reduce levels of greenhouse gases (Ghanbari, 2014).
Changes to land use are a major driver of carbon storage in the terrestrial ecosystem (Chuai et al., 2013).Cellular automata are often used to simulate spatial changes in land use.A cellular automaton is a system that contains a finite number of cells that are located next to each other in a regular and continuous grid which takes the form of a raster (Kwadwo Nti, 2013).When used to model land use, the cells are classified into multiple grid classes according to land use type.This classification includes sorting the cells in groups of cells based on the relationships that exist between them (Sokal, 1974;Anderson et al., 1976;Gregorio and Jansen, 2005).Changes in land use/cover are caused either by natural factors or human intervention.Although human intervention is a long time historical factor in changes to land use/cover, the current rate of change is much higher now than in the past, and this is causing many changes in the environment and eco- logical processes at the local, regional and international levels.Monitoring and modeling regional land use helps to predict future consequences of changes in land use/cover, so that resources can be adequately protected ( Şatır and Berberoglu, 2012).
This study investigates if and/or how changes in land use/cover affect carbon sequestration, which we will regard as an ecosystem service, over time.Furthermore, this study will ascertain if cellular automata are an appropriate simulation method for predicting the levels of carbon sequestration.Specifically, we use cellular automata to simulate different classes of land use/cover.

Study area
The Iranian city of Karaj, located in the region 35.42 • N, 50.50 • E to 35.53 • N, 51. 03 • E (Fig. 1), was selected as the study area.
Karaj was recorded as having a population of around 1.97 million in its 2016 census, making it the fourth-largest city in Iran, after Tehran, Mashhad and Isfahan; in the period from 2006 to 2011 Karaj had a rapid annual growth rate of 3.6 % (SPSDKC, 2012).The city of Karaj became the capital of Alborz Province in 2008 and it has faced rapid development, especially in its urban structures, such that it has been transformed from an agricultural-based city to a city with rapid population growth.High quality green spaces, for example agricultural land and gardens, have rapidly been transformed into industrial, residential and urban service areas.Figure 2 shows land use/cover classes for the years 1985 and 2014.

Simulation of land use/cover changes
Satellite images were obtained from Landsat TM and ETM+ sensors during the period from 1985 to 2014, and land use/cover changes were simulated using the Dinamica EGO 2.4.1 software.Prior to use in the simulations, all data were preprocessed in other applications such as ArcGIS 10.2, ENVI 4.7 and IDRISI Taiga, and all images were checked (number of rows and columns, geographic coordinate system and overlapping border areas), reviewed and modified as needed (Makhdoom et al., 2001).For the cellular automata simulation, the main simulation parameter was the probability of a change in class for each cell.Probability values were calculated using the Markov chain method (Soaresfilho et al., 2012).The probability of a cell's class changing was then used as the rule for a cellular automata, giving a probability map that shows the most probable status for land use/cover classes after a specific period of time.Markov chain probabilities can be improved by adding other auxiliary relations into the cellular automata rule.These auxiliary relations come from regression analysis between land use/cover changes and arbitrary variables such as depth of soil, slope, aspect and distance to rivers and water bodies.
In this study we considered five land use/cover classes: (i) residential, (ii) agricultural, (iii) rangeland, (iv) forest and (v) barren areas.In addition, following auxiliary variables were used: (i) slope, (ii) aspect, (iii) elevation, (iv) distance to rivers and other water bodies, (v) distance to roads, (vi) depth of soil, (vii) salinity of soil, (viii) soil texture and (ix) distance to each other land use/cover class.Before running the simulation based on the proposed model, all the proposed variables were examined for effectiveness on the model; all ineffective variables were removed from the model separately for each transition.The transition matrix was then determined using the Markov chain model.The cellular automata model, which took the initial land use/cover maps as its input, and which used transition matrix as its main rule, was used to simulate land use/cover over a specific period of time.The initial simulation, which was used to validate the model, was performed for the period from 1985 to 2014.A simulation was then performed for the period from 2014 to 2029 to predict the future status of land use/cover classes (Soares-filho et al., 2009;Kwadwo Nti, 2013).
The Dinamica EGO software considers many types of embedded models and it enabled us to introduce data to our model, including land use/cover maps and auxiliary variables.We executed the model for a specified time period and stored the results as maps.The Dinamica EGO software incorporates a number of map comparison tools; we used the "fuzzy reciprocal similarity tool", which locates a window of odd dimensions at the same coordinates on two different maps.Proportional to the relative window sizes, a minimum similarity was calculated using the fuzzy reciprocal similarity method.Briefly, this method uses an enhanced kappa accuracy method which is suitable for circumstances where "small patches" cannot be "moved around" the area under investigation.In our study, we used this tool to compare the predicted and observed maps in order to validate the proposed model.The size of the moving window was varied from 1 to 33 (arbitrary units).The square window, with unity area, had a side length of 30 m and a land area that is equivalent to a square window with an area of 33 and a side length of over 990 m on Earth (Soares-filho et al., 2009(Soares-filho et al., , 2012;;Kwadwo Nti, 2013).

Sampling and calculation of carbon sequestration index per land use/cover class
A total of 20 sampling points were randomly selected, and three vegetation (plant) and soil samples were taken around each point.The sampling points are shown in Fig. 2 (2014 plate).For plants, carbon sequestration was determined using the dry plant weight method, and for soil carbon sequestration was determined using Walkley and Black's method  (Allison, 1965).Plant sampling was usually performed when plants were in their dormant growth season; thus, the plant species that was dominant in the sampling area determined the time of year at which the sample was taken (Karimian, 2009).For samples that were dominated by grass and shrub cover, above-ground and below-ground biomasses were sampled separately.Direct measurements were taken of aboveground biomass by cutting samples and drying them in the laboratory.The dried samples were weighed and the plant basis weight per unit area calculated; weight of above-ground biomass and total biomass per hectare were thereby estimated (Mesdaghi, 2001).The sampling points were made up of plots with specific dimensions (1 m 2 ).Plant litter from each plot was collected as part of the above-ground biomass.Plant litter samples were transported to the laboratory and, after being fully dried, the dry weight of each sample was recorded and the total mass of plant litter per square meter was determined.The below-ground biomass of the sampled species was also sam- pled, usually in 30 × 30 × 30 cm samples.For each sampling point, roots were collected, dried and weighed, so that the root mass per unit area could be calculated.In forest cover, only trees with trunks of at least 10 cm diameter were selected for sampling in order to avoid damaging very young trees.In order to determine the biomass of standing trees, plots of different sizes were used depending on forest density and homogeneity.To do this, plots were randomly selected within forest areas and the number of each tree species was counted.We measured the tree diameter at 1.5 m height in order to determine the tree volume.To determine the bulk tree density, tree branches with specific dimensions were cut, and their dry weight and volume was measured; bulk tree density was determined by dividing weight by volume.Then, by multiplying the volume of trees with the calculated bulk density, biomass weight per unit area was determined.The method was used for sampling forest litter was similar to the method used for herbaceous plants and shrubs (Karimian, 2009).
To study soil characteristics in areas of grass and shrub cover, profile sampling points were drilled in the areas below the plant canopy.Sampling was carried out in the first 30 cm of soil because that is the depth at which microorganisms are  most prevalent and where the majority of the plant roots that were growing were formed.Samples of approximately 1 kg weight were sent to the laboratory for analysis.Sampling in forest areas was also usually carried out in the first 30 cm of soil, in 30 × 30 × 30 cm (or until bedrock was reached) samples; soil was sampled both under trees and in the spaces between trees (Mahdavi et al., 2007).In other land use/cover classes, such as barren areas, only soil samples were taken.Samples of shoots, below-ground plants and plant litter were dried in the open air for several weeks and then their dry weight were determined (Karimian, 2009).Wet samples were dried in an oven at 72 • C for 48 h and the dry weight of each sample was calculated (Gholami, 2008).The weight of organic carbon stored in samples was determined using the combustion method (McDicken, 1997).Samples were ground to a powder and 2 g of each sample was placed in an electric furnace at 550 • C for 5 h (Ghanbari, 2014).After weighing the combustion ashes, the difference between initial weight and weight of the ashes revealed the amount of organic matter in the sample.The weight of organic carbon is calculated using Eq. ( 1), where OC is the weight of organic carbon and OM is the weight of organic matter (Birdsey et al., 2000): Soil samples were dried in the open air, ground and sieved using a 2 mm mesh.The amount of organic matter and organic carbon in the soil samples was determined according to the soil analysis methods of Walkley and Black (Allison, 1965).Specifically, the amount of soil organic matter was determined by measuring the amount of oxide-able organic matter in the soil; after determining the percentage of organic carbon, the amount of soil organic matter and soil organic carbon was calculated by using the Walkley and Black method, and an Erlenmeyer flask (Haghighi, 2003).The percentage of soil organic matter is calculated by Eq. ( 2): where V 1 is the volume of ferrous ammonium sulfate used for the control sample (in mL), V 2 is the volume of ferrous ammonium sulfate used for the sample (in mL), N is the normality of ferrous ammonium sulfate, S is the weight of dried soil sample in the open air and OM is percentage of organic material in the soil sample.Equation (3) gives the percentage of organic matter expressed as the percentage of organic carbon (Haghighi, 2003): where OC is the organic carbon content.Finally, Eq. ( 4) gives the organic carbon content of the soil, in kilograms per hectare (Haghighi, 2003): where Bd is the soil sample density in gram per cubic centimeter, and E is the sampling depth in cm.

Modeling of carbon sequestration changes
Considering the measured data about the amount of carbon sequestration for each land use/cover class in the study www.solid-earth.net/9/735/2018/Solid Earth, 9, 735-744, 2018 area and combining it with the results of the land use/cover change simulation, the rate of carbon sequestration and its variations in accordance with the land use/cover changes was simulated in the future.Thus, if the sampling region has n land use/cover classes, and each class has C i cells, each class has CS i carbon sequestration; if the area of each cell on the ground is S (in hectares), as determined by the cell size of the satellite image, then the total amount of carbon sequestration in the study area can be calculated using Eq. ( 5): (5) In the above equation for each i = 1, 2, . .., n, part of (CS i × C i × S) gives the amount of carbon sequestration in the ith class, in hectares.The amount of carbon sequestration of each class includes amount of carbon sequestration of vegetation shoots, underground organs, all of the vegetation cover, soil and the average for the land use/cover class.3 Results and discussion

Simulation of land use/cover changes
By using the Dinamica EGO software, a cellular model was run for a 29-year period, starting from 1985, to obtain a simulated result for the year 2014.In order to simulate changes in land use/cover classes, the transition matrix between two land use/cover maps between the years 1985 and 2014 was calculated, based on the Markov chain method.By introducing auxiliary variables to the model, the correlation between them and land use/cover classes was calculated using the regression method; one of two correlated variables, or the variables that did not have any impact on the transition matrix, were excluded from the model.Table 1 details the final effective variables.The final probability matrix was then recalculated (Table 2).
As explained in the material and methods section, minimum similarity between observed and simulated result was determined using the fuzzy reciprocal similarity method for the same area regions in order to validate the model (Table 3 and Fig. 3).
The results of simulated land use/cover classes for changes of land use/cover in the period from 1985 to 2014 are shown in Table 4.
Finally, land use/cover classes were simulated for a period of 15 years, from 2014 to 2029.Selected simulation results are presented in Fig. 4 and Table 5.

Sampling and carbon sequestration index calculations
As mentioned in the methods section, 20 sampling points were selected randomly in accordance with literature reports and expert knowledge; vegetation and soil samples were taken at each point in order to measure the level of carbon sequestration.Using the method of Walkley and Black, the amount of sequestered carbon was calculated to provide an index of soil carbon sequestration of the soil (see Table 6 for results).Following the soil experiments, carbon sequestrations of plant samples were measured and, according to the average volume of plants in each plot, the levels of carbon sequestration of plants per plot (and per hectare) were calculated.After the calculation for each sample in same the class by using plant dry weight method, mean values were calculated for each land use/cover class; the mean values were introduced as an index of the average carbon sequestration in vegetation of that class (Table 7).There were different types of vegetation within the overall study area.For example, in agricultural areas most of plants were alfalfa and potatoes.In rangeland areas, most samples included sagebrush, and in the forest areas, plane tree and cedar were the most commonly sampled.
The mean soil and vegetation index values for each land use/cover class (Tables 6 and 7) were combined to give a total index for each class (Table 8).

Simulation of carbon sequestration changes
Carbon sequestration changes based on the changes of the land use/cover classes in the study area were predicted by inserting the carbon sequestration indices and results of land use/cover simulation into Eq.( 5) (Table 9).

Discussion
By measuring the direction of change of land use/cover, we are able to better understand a range of important phenomena and make important predictions into the future (Webster, 2002); therefore this work is invaluable for future environmental planning.In this study, which was based in the city of Karaj, Iran, we determined, in the laboratory, two separate indices of carbon sequestration for soil and vegetation for five land use/cover classes.We simulated the changes in those classes by using satellite images of Landsat TM and ETM+ sensors during the period from 1985 to 2014, and by using the method of cellular automata.Our simulations determined that, in the study period (1985 to 2014), residential and barren areas underwent an annual size increase of approximately 2.63 and 3.76 %, respectively, and that agricultural, forest and rangeland areas decreased annually by approximately 2.37, 9.08 and 2.04 %, respectively.By comparing the results of observed data from 1985 to 2014 with simulated data from 2014 to 2029 we predict that the rate of land use change in the period from 2014 to 2029 will be slower than that observed in the period from 1985 to 2014.According to our simulations, residential space will increase in area by approximately 10 %, agricultural land will reduce by 5 %, forests will increase slightly by approximately 0.3 %, barren areas will reduce by approximately 6 % and rangeland will remain almost constant, with only slight fluctuations in area.Finally, the results suggest that the city of Karaj will experience a reduction in its levels of carbon sequestration due to the impingement of residential areas, which do not contribute to carbon sequestration, into other land use classes.In this study we used a reciprocal fuzzy method to validate the proposed model by comparing observed and simulated land use/cover maps from 2014.Our results show that observed and simulated results will converge over a 15 year period as the size of the compared areas increases.More specifically, simulated results will converge with observed data for regions with an area of at least 900 m 2 with a probability of > 30 %, and for regions with an area of 26 100 m 2 with a probability of > 70 %.Sheng et al. (2012), Soares-fihlo et al. (2012, 2013) and Kwadwo Nti (2013) also used cellular automata and Markov chain with fuzzy analysis methods to simulate and validate predicted changes in land use/cover use; this further validates the results of the present study.
The carbon sequestration levels were determined for both vegetation and soil samples at each sampling point.Carbon sequestration in soil samples taken from agricultural land was much higher than equivalent samples taken from forest and rangelands.Although there is no current single data source that details the impact of different land management practices on carbon sequestration in agricultural land (World Bank, 2012), it has been demonstrated that carbon sequestration changes according to plant growth season (Forge, 2001); www.solid-earth.net/9/735/2018/Solid Earth, 9, 735-744, 2018 monitoring of carbon sequestration, therefore, remains important.In this study, soil samples were collected at the start of the growing season, in April; this is a season during which farmers use many types of fertilizers, both natural and chemical.Carbon sequestration levels were found to be high in the agricultural classes, especially directly after the application of fertilizers.Indeed, the amount of fertilizer present in the soil was measured so that, if desired, its effect could be removed from the data by using only vegetation-derived data; we decided not to remove those data so that we were able to show the overall impact of land use/cover change simulation on carbon sequestration.The results of other studies confirm that the model used in this study is well suited to simulate land use changes (Al-ghamdi, 2012;Guan and Clarke, 2010;Jantz et al., 2010;Kwadwo Nti, 2013;White andEngelen, 1993, 1997;Deadman et al., 1993;Itami, 1994;White, 1998;Li and Yeh, 2002;Barredo et al., 2003Barredo et al., , 2004;;Al-Ahmadi et al., 2008).Furthermore, the results of this study agree with the findings of Soares-filho et al. ( 2012), who used cellular automata and Markov chain tools with Dinamica EGO software to study carbon sequestration in agricultural land and forests in Brazil.
It is important to note that the results of this study are specific to urban areas, and that they can, therefore, be applied to urban growth over time.For our study area, as detailed in Tables 4 and 5, some land use/cover classes, such as agricultural, decrease in size due to increases in residential requirements.Thus, our results cannot be extrapolated to all landscape types beyond urban.Furthermore, it is important to note that the results of our study indicate that the use of fertilizers resulted in levels of carbon sequestration that were elevated in agricultural areas compared to forest areas.We decided to include data from soil that had been fertilized because (i) the traditional fertilization methods used by the farmers meant that the amount of fertilizer used was not recorded and (ii) agricultural land forms a significant part of carbon sequestration capacity and, therefore, should not be excluded from the study.

Conclusions
Our results indicate a number of opportunities to include carbon sequestration in simulations of land use change in fast growing cities; the methods used in this study can be applied elsewhere as an "ecosystem service" for urban planning based on land use changes.Based on our results, which demonstrate rapid changes in residential areas, we recommend that city managers take effective measures to increase residential green spaces which are planted with appropriate plant species.We also recommend that comparative studies are performed to confirm that our methods can be used elsewhere; these studies should ascertain whether the results are affected by variables that were not possible to isolate in the present study, such as population density, crop yields, or biological and physical soil characteristics.
Data availability.All used data in this research were provided by the authors.These data are divided into two parts: the first part consists of the amount of carbon sequestration for sampling points that were measured in the soil laboratory at the faculty of natural resources at the University of Tehran, and the second part consists of satellite image processing that was done in the laboratory of environment at the faculty of natural resources at the University of Tehran.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.The location of the city of Karaj in Iran.

Figure 3 .
Figure 3. Maximum and minimum similarities between simulated and observed maps for the year 2014.

Figure 4 .
Figure 4. Predicted land use/cover classes in different years.

Table 1 .
Summary of final effective variables.

Table 3 .
Minimum similarities between observed and simulated results for the year 2014.

Table 4 .
Changes of land use/cover in the period from 1985 to 2014 (in hectares).

Table 5 .
Predicted changes to land use/cover in the period from 2019 to 2029 (in hectare).

Table 6 .
Carbon sequestration of soil by land use/cover class.

Table 7 .
Carbon sequestration in vegetation per land use/cover class, using dry plant weight method.

Table 8 .
Total carbon sequestration per land use/cover classes.

Table 9 .
Simulated carbon sequestration of the study area from 2014 to 2029 (millions of tons).