The Carbon Dynamics of Dry Tropical Afromontane Forest Ecosystems in the Amhara Region of Ethiopia

Forest degradation due to land use change is a severe problem in Ethiopian Afromontane Forests. We investigated such degradation effects by comparing degraded agricultural land (previously covered with forest) with neighboring natural forests, 40 to 50 years after conversion. We selected four different study areas to cover the eco-climatic conditions of the Amhara region in Northwestern Ethiopia. For a paired-stand comparison we collected soil data on both land use types. We calculated forest biomass to evaluate the biogeochemical-mechanistic ecosystem model Biome-BGC, which is used as a diagnostic tool to assess the site and management impacts on productivity as well as ecosystem carbon and nitrogen accumulation. We applied Biome-BGC to assess rehabilitation options on such degraded land. Afromontane forests in the highlands of Ethiopia showed high soil C stocks, resulting from long lasting biomass accumulation. Removing the tree cover and converting forest areas to cropor grassland, has led to a loss of 40–85% of the soil C stocks and thus a loss in soil fertility within only 40 to 50 years. Rehabilitation efforts by replanting trees will improve soil fertility, but may require over a thousand years to achieve a similar level in biomass and soil fertility versus the situation prior to the land use change.


Introduction
Tropical forests account for 50% of the Earth's total plant biomass [1] and store 40% (428 Gt C) of the terrestrial carbon [2]. Tropical forests are threatened by deforestation due to land use changes. A good example is Ethiopia, which experienced a decline in forest covered land area from 40% to less than 10% during the last century [3]. The removal of the forest cover has resulted in severe soil erosion, species extinction, as well as reduction in productivity and carbon sequestration [4]. The main reasons for the rapid land use change is the expansion of agricultural land associated with the demand for fuelwood for a fast growing population [5]. Ethiopia has a population of about 100 million people and an estimated annual increase of 3% or 3 million people, 73% of them working in the agricultural sector [6]. For the rural population, it is increasingly difficult to cover the supply of fuel wood for heating homes and cooking, since only small forest areas are left, and the remaining forests are largely protected. In most of the cases, these forests remained because they are located next to monasteries and are preserved as "church forests". However, these "church forests" are not untouched, since we can observe strong evidence for forest management impacts.

2.
Assess the degradation effect of land use change from forest to agricultural land. 3.
Analyze the rehabilitation processes following afforestation on degraded agricultural land using Biome-BGC.

Study Areas
We selected four study areas covering different agro-ecological zones within the Amhara region in Northwestern Ethiopia (Figure 1). The four sites were selected to represent different forest ecosystems ranging from low (860 m a.s.l.) to high (2500 m a.s.l.) altitudes, a mean total annual precipitation from 1080 mm to 1708 mm, and a mean annual temperature from 17.1 • C to 26.1 • C ( Figure 2). Each study area covers two different land use forms located next to each other: (i) natural forests and (ii) agricultural land. The agricultural land was established on former natural forest about 40 to 50 years ago [32]. This paired-plot approach enabled us to directly address the impact of land use changes, as climatic and site conditions can be assumed having been similar.
The four study areas are located in Gelawdiwos, Injibara/Katassi, Mahibereselassie and Taragedam. The forest areas in Gelawdiwos, Taragedam and Mahibereselassie are owned by the Ethiopian Orthodox Tewahedo Church, whereas the forest in Katassi is owned by the local government. Gelawdiwos represents the growing conditions in the highlands of the Amhara region. The forest composition is mostly dominated by evergreen highland tree and shrub species, such as Chionanthus mildbraedii (Gilg & G.Schellenb.) Stearn, Euphorbia abyssinica J.F.Gmel. and Ekebergia capensis Sparrm. Katassi and Taragedam are both located in mid-high altitudes. The Katassi forest is dominated by the mid-altitude Afromontane tree and shrub species Albizia schimperiana Oliv., Prunus africana (Hook.f.) Kalkman and Croton macrostachyus Hochst. ex Delile. Similarly, Taragedam is composed of many Afromontane forest species, including Juniperus procera Hochst. ex Endl., Teclea nobilis Delile and Schefflera sp. J.R.Forst. & G.Forst.; Mahibereselassie is the typical ecosystem in the lowlands of the Amhara region and contains deciduous and semi-deciduous woodland and savanna grassland tree and shrub species such as Sterculia setigera Delile., Boswellia papyrifera Hochst. and Terminalia laxiflora Engl.      (Walter and Lieth, 1960) for the four study areas (dotted areas indicate dry periods, hatched area humid periods and black areas wet periods). The graphs also show elevation, climate data recorded period, the mean annual temperature and the mean total annual precipitation.

Forest Data
Forest data in the four study areas were collected based on a systematic grid sampling approach ( Figure 1, Table 1). The grid size of the systematic sampling within each of the four study areas differed according to the size of the forest. Slope, aspect, latitude and elevation were recorded for each plot.
Since average stand density varied by study area but forest inventories require a certain average minimum number of trees per plot, we varied the plot size by study area. In Gelawdiwos and Taragedam, information for trees with a DBH (diameter at breast height) >10 cm was collected on plots with a 10 m radius, while for Mahibereselassie and Katassi a plot size radius of 15 m was chosen. On these plots, the following tree information was recorded: tree species, tree location, DBH, crown width in the four cardinal directions, tree height and height to live crown base. In addition, increment core data were collected by major species groups on each sampling plot and for each "central stem". The term "central stem" refers to the 60th percentile of the DBH distribution and follows a suggested procedure by [36], which showed that the central stem represents the "mean tree" (=quadratic mean DBH) on a given plot. Walter/Lieth climate diagram (Walter and Lieth, 1960) for the four study areas (dotted areas indicate dry periods, hatched area humid periods and black areas wet periods). The graphs also show elevation, climate data recorded period, the mean annual temperature and the mean total annual precipitation.

Forest Data
Forest data in the four study areas were collected based on a systematic grid sampling approach ( Figure 1, Table 1). The grid size of the systematic sampling within each of the four study areas differed according to the size of the forest. Slope, aspect, latitude and elevation were recorded for each plot. Since average stand density varied by study area but forest inventories require a certain average minimum number of trees per plot, we varied the plot size by study area. In Gelawdiwos and Taragedam, information for trees with a DBH (diameter at breast height) >10 cm was collected on plots with a 10 m radius, while for Mahibereselassie and Katassi a plot size radius of 15 m was chosen. On these plots, the following tree information was recorded: tree species, tree location, DBH, crown width in the four cardinal directions, tree height and height to live crown base. In addition, increment core data were collected by major species groups on each sampling plot and for each "central stem". The term "central stem" refers to the 60th percentile of the DBH distribution and follows a suggested procedure by [36], which showed that the central stem represents the "mean tree" (=quadratic mean DBH) on a given plot.
In addition, for trees with DBH <10 cm and a height >1.5 m, the same set of tree variables (except the increment cores) was collected on subplots with a radius of 4 m. All trees (from plots and subplots) with a DBH ≥2.5 cm were used for deriving vegetation carbon and carbon increment. Missing tree heights H were derived from the Petterson height-diameter curve [37]: where the parameters a and b were derived for the twenty most important tree species in the study areas [38]. Above ground carbon vegetation C veg was calculated according to the following formulation [39]: where DBH and H are as previously defined and is the mean wood density of African tree species set equal to 0.614 g cm −3 [40]. The carbon content was assumed to be 50% of the biomass [41][42][43]. The annual NPP (g C m −2 year −1 ) of the woody part of each tree was calculated as the mean annual woody carbon increment for the 10 year period 2005-2014. The increment was derived from the increment cores of the central stem [38] since the growth rate of the central stem is proportional to the growth rate per hectare, assuming that the stem number remains equal within the growth period. The total annual terrestrial NPP for each inventory plot is the sum of woody carbon increment, the litterfall and the root increment (fine root production and coarse root increment) [44]. Since no litterfall and root data are available, we followed [45], who suggest a wood increment of about 40% of the total carbon increment. The remaining 60% are the litterfall and the root increment (30% each). For further details regarding the sampling method and details related to the calculation of the standing timber volume and volume increment rates, we refer to [38]. Summary statistics of the forest data by study area are shown in Table 1. Note that the number of plots in this study are slightly lower versus [38] because we eliminated plots with no tree coverage at the edge of the sampled forests.

Forest Area
Soil samples from the forested areas were systematically taken on the same general grid as the forest biomass samples, at about one fourth of the forest sampling grid points (see Figure 1 and Table 1). The litter layer was removed and soil probes were taken with a soil corer (Pürckhauer type, inner diameter 30 mm). The soil samples were taken in 0-10 cm, 10-20 cm, 20-30 cm, and 30-50 cm soil depth and stored in plastic bags. The probes were sieved to 2 mm and soil organic carbon, nitrogen, and texture were determined in the laboratory in Vienna.
The relative weight proportion of the soil organic matter was calculated as the loss of ignition (LOI), heating the soil samples in a muffle furnace [46]. After oven-drying the samples at 105 • C until a constant weight was achieved, the organic matter was combusted to ash at a temperature of 450 • C for 4 h. The LOI was then calculated using the following equation: where LOI 450 represents LOI at 450 • C (as a percentage), DW 105 represents the dry weight of the sample before combustion and DW 450 the dry weight of the sample after heating to 450 • C. The weight loss is proportional to the amount of organic carbon contained in the sample. Carbon stocks were determined on a weight to area basis (kg C m −2 ). The carbon concentration was corrected for bulk density and soil volume. The total carbon TC (kg C m −2 ) was calculated as: where C is carbon concentration (%), BD is bulk soil density and V is volume of soil for each sample plot. Bulk soil density was determined from samples taken separately on the side of each soil carbon sampling point. Soil texture was determined by a combination of wet sieving and sedimentation methods [47]. Texture classes were defined according to the USDA particle size classes sand (2.0-0.053 mm), silt (0.053-0.002 mm), and clay (<0.002 mm). Soil depth as needed for Biome-BGC was obtained from the regional soil depth map provided by Amhara Design and Supervision Works Enterprise (ADSWE) [48]. Based on the ADSWE soil depth ranges and the local experience, we defined for the shallow lowland site (Mahibereselassie) a range of 40 cm and the deeper highland sites a range of 1 m ( Table 2) to cover uncertainty from natural variation in soil depth.

Agricultural Area
In addition to the forested areas, we obtained soil data from neighboring agricultural sites (cropland and grazing land). A total of 40 sample plots, 10 plots for each study area, were taken on transects at 50 m/100 m intervals. All soil samples were taken at the same time as the forest soil samples. For further details, we refer to [32].

Atmospheric Data
We obtained downscaled 1 km × 1 km daily meteorological data, i.e., daily minimum and maximum temperatures ( • C), daily precipitation (mm) from Sisay et al. [49]. The data cover a 32 years period ranging from 1979 to 2010 and have been downscaled from global precipitation and temperature data on a 1 km × 1 km resolution for the whole Amhara region. The downscaled precipitation was further improved based on the difference in mean annual precipitation between the downscaled dataset and the closest climate station to each study area. Each daily value of the downscaled precipitation was multiplied with the ratio of observed mean annual precipitation to downscaled mean annual precipitation, to ensure correct precipitation sums. Solar radiation (Ws −1 ), vapor pressure deficit (VPD, Pa) and day-length (from sunrise to sunset) were generated from daily temperature and precipitation using algorithms implemented in the Mountain Climate Simulator (MT-CLIM) [50,51].
Current and preindustrial atmospheric CO 2 concentration and nitrogen deposition, as they are needed for the BGC simulations, were taken from two global datasets published by [52,53]. As current nitrogen deposition values, we used the two different values given by [53] for the entire study area and their mean ( Table 2) to address the high variation and uncertainty stemming from the coarse resolution of the dataset. The preindustrial value of atmospheric CO 2 concentration was set to 278 ppm [52,54] and to 0.129 g N m −2 year −1 for Nitrogen deposition [53]. Table 2. Soil depth classes and atmospheric nitrogen deposition levels for the evaluation of Biome-BGC. For each simulation, we used three different soil depths (compare Section 2.3.1) and 3 different nitrogen deposition rates (min/mean/max of values from [53]) to cover the potential variation. The soil depth numbers are taken from the regional soil depth map provided by Amhara Design and Supervision Works Enterprise (ADSWE) [48] and the nitrogen deposition levels come from [53].

The Ecosystem Model
For the simulation of forest carbon storage (objectives 1 and 3), we employed the biogeochemical-mechanistic ecosystem model Biome-BGC [55,56], which simulates flux dynamics of the water, carbon, and nitrogen cycle for a terrestrial ecosystem on a daily time step. For this study, we used Biome-BGC version 4.1.2 [29,51] with an improved self-initialization algorithm to mimic undisturbed ecosystems [57,58] and forest management routines e.g., historic land use changes, thinning, clear-cut and planting [30]. Biome-BGC is fully prognostic and does not require detailed forest data for model initialization. It is initialized by a so called "spin-up" run [27,29,57,59], which mimics the accumulation of soil and vegetation carbon and nitrogen under constant conditions (available climate data recycled, constant preindustrial CO 2 and nitrogen deposition) until a "steady-state" is reached (typically after few thousands of years). Since CO 2 concentration and nitrogen deposition have increased and most forests have experienced some form of forest management, which may have led to forest degradation, the next step is typically to perform "modern-forest-simulations" addressing these changes to realistically mimic the current ecosystem fluxes for a given forest [27,30].
Key processes within Biome-BGC are the photosynthetic carbon assimilation, autotrophic respiration (maintenance respiration plus growth respiration), allocation among different plant compartments, mortality, litterfall, mineralization processes, nutrient and soil water uptake, rainfall interception, soil water infiltration and storage, runoff and evapotranspiration. Leaf area determines potential forest productivity as well as rainfall partitioning (interception, through fall) and evaporation from soil. Carbon assimilation is calculated with the Farquhar photosynthesis routine [60].
Growth respiration is a fixed fraction of the amount of carbon allocated to different plant structural compartments [61], whereas maintenance respiration depends on the nitrogen content of the tissue [62]. CO 2 availability for photosynthesis is regulated by the stomata, where the CO 2 influx depends on maximum conductivity, and stomatal closure is limited by radiation, moisture of the atmosphere and water content in the leaves, which is controlled by the soil water potential. Soil water potential is calculated from the soil water and texture following [63,64]. Key model outputs are Gross Primary Productivity (GPP), Net Primary Productivity (NPP; GPP minus autotrophic respiration), and carbon storage in the living biomass, the litter, coarse woody debris and the soil.

Eco-Physiological Parameters
Within Biome-BGC, biomes or species are defined by 39 "Eco-physiological parameters", which characterize a forest ecosystem (and are stored in the epc-file). The original Biome-BGC model was developed for seven different biome types (evergreen broadleaf, deciduous broadleaf, evergreen needle, deciduous needle, shrubs, and C 3 grass (temperate or cool-season grasses, using C 3 carbon fixation) and C 4 grass (tropical or warm-season grasses, using C 4 carbon fixation) [65]. Over the years some regional or species-specific parameter sets were published [66,67].
For our study, we selected the parameter set proposed for tropical forests by [66]. For the driest study sites located in the lowland area, Mahibereselassie, we changed the evergreen phenology, because most tree species in the lowland area are dry-deciduous or semi-deciduous and the flushing of leaves as well as litterfall is controlled by the precipitation pattern. Budburst starts at the beginning of the rainy season in April, and the annual litterfall is in December. According to these regional conditions, we changed the date of the vegetation period in the epc-file and set the start date of the growing season in Mahibereselassie to day 90 and the end date (final dropping of the leaves) to the last day of the year (day 365).

Current Forests
For the simulation of the current forests in the four study areas with Biome-BGC, we used the atmospheric, soil (texture, depth) and eco-physiological input data as described in Sections 2.3, 2.4 and 3.2, respectively. The variations in soil depth and nitrogen deposition we employed are shown in Table 2.
We started with the spin-up run followed by a "realistic" land use history scenario ( Figure 3). For the spin-up, the available climate records were used repeatedly, after having been corrected based on the assumption that Ethiopian highlands had experienced long lasting cooling periods during the Holocene [68][69][70]. We addressed this historic cooler climatic periods in our spin-up runs by reducing the daily mean temperature of the available daily climate data by 3 • C and by increasing the daily precipitation by 20% [69]. For our study, we selected the parameter set proposed for tropical forests by [66]. For the driest study sites located in the lowland area, Mahibereselassie, we changed the evergreen phenology, because most tree species in the lowland area are dry-deciduous or semi-deciduous and the flushing of leaves as well as litterfall is controlled by the precipitation pattern. Budburst starts at the beginning of the rainy season in April, and the annual litterfall is in December. According to these regional conditions, we changed the date of the vegetation period in the epc-file and set the start date of the growing season in Mahibereselassie to day 90 and the end date (final dropping of the leaves) to the last day of the year (day 365).

Current Forests
For the simulation of the current forests in the four study areas with Biome-BGC, we used the atmospheric, soil (texture, depth) and eco-physiological input data as described in Sections 2.3, 2.4 and 3.2, respectively. The variations in soil depth and nitrogen deposition we employed are shown in Table 2.
We started with the spin-up run followed by a "realistic" land use history scenario ( Figure 3). For the spin-up, the available climate records were used repeatedly, after having been corrected based on the assumption that Ethiopian highlands had experienced long lasting cooling periods during the Holocene [68][69][70]. We addressed this historic cooler climatic periods in our spin-up runs by reducing the daily mean temperature of the available daily climate data by 3 °C and by increasing the daily precipitation by 20% [69]. Next, the carbon and nitrogen stocks resulting from this self-initialization had to be adjusted for forest management impacts [28]. Note that even the simulated protected "church forests" cannot be considered as pristine forests because some management impacts (e.g., harvesting, etc.) were evident. For all four forest areas, cutting of selected high valuable trees and the removal of  Next, the carbon and nitrogen stocks resulting from this self-initialization had to be adjusted for forest management impacts [28]. Note that even the simulated protected "church forests" cannot be considered as pristine forests because some management impacts (e.g., harvesting, etc.) were evident.
For all four forest areas, cutting of selected high valuable trees and the removal of suppressed trees for fuelwood was reported (from interviews, documents, [32]). Based on these indications of intensive and uncontrolled extractions, we assumed an annual harvest of 5% of the biomass, as harvest must be assumed to be clearly more than the 3% mortality from the eco-physiological parameter set [66]. As related to the annual harvest, the same share of leaves and fine roots was transformed into litter and the same share of coarse roots was moved to the woody debris carbon pool.

Rehabilitation Pathways
We also used Biome-BGC to assess the rehabilitation potential of afforestation on degraded agricultural land and the time it may take until a new "equilibrium" or "steady state" in the accumulated biomass and soil carbon stocks is achieved. After manually setting the soil carbon and nitrogen values in the model equal to the recorded values on the degraded agricultural land (see Table 3) and assuming a "standard BGC planting" of forests (using the previous tropical forests epc-file), the simulations were started. For simplicity, the current climate data were repeatedly used without assuming any future changes in climate (because reasonable climate change projections for the next few thousand years do not exist). Since undisturbed forest development seems to be unrealistic due to the enormous demand for fuelwood, we employed three different management scenarios: (i) no management, (ii) 2.5%, and (iii) 5% annual extraction of timber. The corresponding amount of leaves and fine roots was transferred to litter pool, and coarse roots to coarse woody debris, respectively. Table 3. Summary of the available soil data; The grid data refer to the soil carbon, nitrogen and Sand/Silt/Clay information collected on parts of the forest inventory grid points; The transect data are based on [32] and refer to transects with 10 samples taken every 50-100 m, placed in the agricultural area of our study areas.

Evaluation of Biome-BGC for Ethiopian Afromontane Forests
Observed vegetation carbon stocks, soil carbon and NPP differed markedly among the four study areas. Gelawdiwos, the study area highest in elevation showed the highest productivity and carbon storage, whereas Mahibereselassie, located in the dry lowlands showed the lowest productivity and carbon storage. Taragedam with the second lowest precipitation also showed the second lowest productivity and carbon storage.
Biome-BGC simulations followed this general pattern. Note that the assumed management with 5% extraction annually was identical for the performed simulations in the four study areas. Figure 4 shows the median, the 25th and 75th percentile values of predicted versus observed vegetation and soil carbon. The observed vegetation carbon of the four study locations ranged from 1 t C ha −1 to 187 t C ha −1 . The corresponding Biome-BGC predictions of the carbon content of the vegetation carbon for nine combinations of soil depth and nitrogen deposition levels ranged from 11 t C ha −1 to 67 t C ha −1 . The observed soil carbon ranged from 14 t C ha −1 to 253 t C ha −1 and the modeled numbers ranged between 16 t C ha −1 to 242 t C ha −1 .

Land Use Change Effects
We compared the soil carbon and nitrogen content for the forested areas and the neighboring agricultural land (Table 3). For three sites (Gelawdiwos, Katassi and Mahibereselassie), the available soil information came from crop land, while in Taragedam, the soil data came from grazing land since no crop land was available in the area [32]. All sites showed significantly lower soil carbon and nitrogen stocks on the agricultural land compared to the forest, except the dry lowland study area Mahibereselassie. The Mahiberseselassie site is not a typical forest area, since trees do not form closed forest stands due to limitations in precipitation. In addition, the "forest area" experienced intensive grazing and regular grass burning impacts [32]. The most dramatic reduction in soil carbon stocks was observed in Katassi, where 85% of the soil carbon was lost. These differences are the result of approximately 40 to 50 years of intensive farming since the conversion of forest to agricultural land, including the typical erosion effects during the rainy seasons in the highlands of Ethiopia.

Rehabilitation
Next, we evaluated the soil amelioration potential of afforestation of these degraded soils by simulating afforestation of the three severely degraded sites with Biome-BGC. After degradation, it could take a few decades to over 2000 years until a new steady state in nitrogen and carbon pools would be reached ( Figure 5). The steady state of soil carbon for the three highland study areas lie between 100 t C ha −1 and 125 t C ha −1 for 5% annual extraction, between 130 t C ha −1 and 155 t C ha −1 for 2.5% annual extraction and 200 t C ha −1 to 235 t C ha −1 for no management. For our simulation, we used the current climate conditions, since employing climate change scenarios for such a long projection period is unrealistic.
During recovery, two marked changes in accumulation velocity of vegetation and soil carbon can be observed. The first marked change in carbon accumulation velocity takes place after some decades, when the system reaches its equilibrium between the annual biomass accumulation and the biomass loss through mortality and harvesting. Net soil and vegetation carbon build up continue on a slower rate past this point, because the system is still accumulating nutrients (i.e., nitrogen from nitrogen deposition and nitrogen fixation). Only much later (some hundreds or thousands of years NPP derived from the recorded data, ranged between 4 g C m −2 year −1 to 1567 g C m −2 year −1 (Figure 4). Corresponding Biome-BGC predictions were within a range of 186 g C m 2 year 1 to 1095 g C m −2 year −1 . Except for Mahibereselassie, predicted NPP was in the range of observed NPP.

Land Use Change Effects
We compared the soil carbon and nitrogen content for the forested areas and the neighboring agricultural land (Table 3). For three sites (Gelawdiwos, Katassi and Mahibereselassie), the available soil information came from crop land, while in Taragedam, the soil data came from grazing land since no crop land was available in the area [32]. All sites showed significantly lower soil carbon and nitrogen stocks on the agricultural land compared to the forest, except the dry lowland study area Mahibereselassie. The Mahiberseselassie site is not a typical forest area, since trees do not form closed forest stands due to limitations in precipitation. In addition, the "forest area" experienced intensive grazing and regular grass burning impacts [32]. The most dramatic reduction in soil carbon stocks was observed in Katassi, where 85% of the soil carbon was lost. These differences are the result of approximately 40 to 50 years of intensive farming since the conversion of forest to agricultural land, including the typical erosion effects during the rainy seasons in the highlands of Ethiopia.

Rehabilitation
Next, we evaluated the soil amelioration potential of afforestation of these degraded soils by simulating afforestation of the three severely degraded sites with Biome-BGC. After degradation, it could take a few decades to over 2000 years until a new steady state in nitrogen and carbon pools would be reached ( Figure 5). The steady state of soil carbon for the three highland study areas lie between 100 t C ha −1 and 125 t C ha −1 for 5% annual extraction, between 130 t C ha −1 and 155 t C ha −1 for 2.5% annual extraction and 200 t C ha −1 to 235 t C ha −1 for no management. For our simulation, we used the current climate conditions, since employing climate change scenarios for such a long projection period is unrealistic.
1 Figure 5. Soil carbon (first row) and vegetation carbon trend (second row) attained on degraded sites after afforestation, without human intervention (left) and with human intervention (2.5% thinning, center and 5% annual thinning, right).
During recovery, two marked changes in accumulation velocity of vegetation and soil carbon can be observed. The first marked change in carbon accumulation velocity takes place after some decades, when the system reaches its equilibrium between the annual biomass accumulation and the biomass loss through mortality and harvesting. Net soil and vegetation carbon build up continue on a slower rate past this point, because the system is still accumulating nutrients (i.e., nitrogen from nitrogen deposition and nitrogen fixation). Only much later (some hundreds or thousands of years later), equilibrium between nutrient accumulation and nutrient loss is reached, and soil and vegetation carbon stocks stay more or less stable.

Discussion
Biogeochemical-mechanistic modeling is a powerful diagnostic tool to assess the carbon and nitrogen fluxes within our study areas. Although the species composition of the natural Afromontane forests in Ethiopia is very rich and diverse, the eco-physiological parameter settings suggested by [66] are applicable within the area. After adjusting for the length of the growing season (on the site with predominant deciduous vegetation cover) and addressing harvesting impacts, the above ground biomass, as well as soil carbon and NPP predictions derived with Biome-BGC revealed consistent estimates versus field observations (Figures 4 and 5). This suggests that our chosen approach mimics the ecosystem fluxes (i.e., carbon, nitrogen and water) as well as potential harvesting impacts correctly for Afromontane ecosystems of the Amhara region of Northwestern Ethiopia. This is an important finding because forests cover about 251,000 ha or 2% of the total land area (15.7 million. ha) in Amhara region [38], but no systematic forest monitoring or empirically driven forest growth and yield model for ensuring sustainable forest management is in place. Considering the enormous pressure of a fast-growing population on the remaining forest resources, this lack of information is one of the key obstacles for sustainable forest management in the region. The option we tested in this study to circumvent this lack of information is the application of biogeochemical-mechanistic modeling theories.
Although the model results fit with the observed data (Figures 4 and 5), the limitations of the modeling approach versus growth and yield models is the level of abstraction in the general modeling approach. This makes forest stand specific predictions, including single tree mortality or specific harvesting operations, difficult [67]. This may be one of the reasons that the model predicted higher vegetation carbon versus the observations in Taragedam (Figure 4), as the study sites in Taragedam experienced stronger human interventions and cattle grazing versus the other study sites [71,72].
The key problem within the highlands of the Amhara region is land use change (the conversion of forests into agricultural area) to feed a fast growing population [4,73,74]. An enormous loss in soil fertility due to erosion is associated with this reduction in forest area. We observed that 40 to 50 years after converting forests into agricultural land, the carbon and nitrogen content in the soil declined by 40% to 85% of the remaining neighboring forested areas ( Table 3). The dramatic deterioration of the soils within this short period of time is also reported by other studies showing similar trends [75,76]. The underlying parent material in Amhara is granite with an extremely low release of nutrients from weathering. Removal of the forest cover leads to enormous soil erosion during the heavy and intensive rainy season [34], associated with a removal of the top soil layer. Once the forest vegetation and the top soil layer are gone, the key source for nutrients is gone [32]. Future productivity is therefore jeopardized by two effects, i.e., loss in nutrients through vegetation clearing and erosion, and the low rates of replenishing of the lost nutrients from weathering of the rock and from nitrogen deposition and fixation [77].
One option to stabilize and eventually ameliorate productivity is to replant forests on such degraded sites. Any permanent land cover will decrease further erosion and will lead to nitrogen and litter input due to litterfall and consequently to a recovery of the carbon and nutrient stocks, and the water retention capacity. The rate of carbon pool accumulation varied by harvesting intensity ( Figure 5). Regardless of the employed harvesting scenario, soil fertility will increase again, but it will stabilize at different levels according to the harvesting intensity. Remarkably, the soil carbon equilibrium levels for the two management scenarios lie below the current forest soil carbon stocks. This indicates that the current higher forest soil carbon stocks ( Figure 4) are a legacy of a former 'better' climate (compare Section 3.2, climate assumption for the spin-up: 20% higher precipitation and 3 • C lower temperature in comparison to today). Thus, although afforestation is strongly recommended since it has a high ameliorating potential of degraded soils, it is questionable if the mistakes from the past can ever be fully repaired under today's climatic conditions. We did not use climate change scenarios for above mentioned reason. However, a scenario of more irregular and intensive precipitation events [78] would increase erosion and slow down the soil recovery especially in the already water-limited lowland areas of Amhara region. In any case, the importance of avoiding deforestation in the first place, is the key message which needs to be conveyed to any government in tropical developing countries around the globe.

Conclusions
With this study, we provided insights into how land use changes-converting natural forest areas into agricultural land-affect the soil fertility and thus the productivity expectations of Afromontane forests in the Amhara region of Northwestern Ethiopia. Within 50 years, soil fertility expressed by the carbon and nitrogen content in the soil has dropped by half or more versus areas covered with forests; and rehabilitation may take hundreds or thousands of years if soil fertility prior to the land use change were to be reached. Countries such as Ethiopia, where no forest inventory system exists, urgently needs a tool to assess carbon storage in forest ecosystems and analyze land use impacts as well as associated productivity changes. We could show that the ecosystem model Biome-BGC can serve as such a diagnostic tool, owing to the model's mechanistic representation of natural biogeochemical processes based on the interactions of climate, site, the simulated vegetation type and management. The model could, therefore, be utilized for large-scale forest productivity estimates and afforestation planning for selecting the most suitable sites in Amhara, Ethiopia. also go to Amhara Regional Agricultural Research Institute in hosting the research project and providing field cars. We thank Hadera Abraha, Khlot Gebrehana, Sibhatu Abera and Tesfaye Teklehaymanot for helping in the data collection. We thank the editor and the three reviewers which helped to improve the manuscript.
Author Contributions: B.B. analyzed the data, did the scenario analysis and wrote the first draft of the manuscript, E.P. helped with model runs and their interpretation and improved the manuscript, K.S. helped with the data collection, D.A. helped with the soil data collection and laboratory analysis, and H.H. designed the project, coordinated and supervised the data collection, the analysis, and the manuscript writing.

Conflicts of Interest:
The authors declare no conflict of interest.