Land use change in the Ecuadorian páramo: The impact of expanding agriculture on soil carbon storage

ABSTRACT

The Andean páramo is notable for high soil carbon storage and its contribution to ecosystem services. However, the páramo’s ability to maintain high soil carbon levels is threatened by land use change from tussock grassland and shrublands to agricultural uses. A chronosequence study was conducted in the páramo around Quito, Ecuador, to determine the rate of soil carbon loss from traditional fallow agriculture. In parallel, a land use and land cover classification of Landsat images was used to measure the change in agricultural areas between 1991 and 2017. There was a significant negative relationship between the time since initial cultivation of a field and soil C: Older agricultural sites had significantly less C compared to natural ecosystems due to an average loss of 0.045 percent soil C per year. Undisturbed sites had significantly more soil C than cultivated sites but not pastures or fallow fields, indicating that cultivation is the most detrimental stage of the fallow agricultural cycle for soil C storage. There was an 838 percent increase in cultivated land between 1991 and 2017 but a 10 percent decrease in pastures, indicating a trend away from traditional regenerative agriculture toward land use types that lead to substantial losses in soil carbon.

KEYWORDS:

• Andes
• ecosystem services
• fallow agriculture
• Andosols

Introduction

Páramo ecosystems, a diverse set of high elevation vegetation types from the Northern Andes, offer a striking example of the disproportionate role that mountain ecosystems have in the provision of ecosystem services (Lazo et al. 2019). Páramo ecosystems account for less than 1 percent of the total land in Venezuela, Colombia, Peru, and Ecuador, yet millions of people from these countries depend on páramo ecosystems for the provision of water for human consumption and sanitation, and large portions of their economy rely on that same water supply for irrigation and generation of energy (International Union for Conservation of Nature 2002). However, though freshwater is arguably one of the most important benefits from the páramo, the ecosystem’s ability to provide such a service stems from its cool climate and the unique characteristic of its soils.

Páramo soils, which are known for their high C storage capacity, are generally classified into Andosols, Entisols, Inceptisols, and Histosols (Batjes 1996). The cool climate in this environment slows decomposition, leading to the accumulation of large amounts of soil organic matter. High carbon storage capacity is also partly due to the high Al concentrations that form stable Al complexes with soil organic matter and also offer protection from microbial activity stemming from Al toxicity and low pH (Tonneijck et al. 2010). As a result, reported soil C concentrations range between 239 MgC·ha⁻¹ in well-drained upland soils and 2,100 MgC·ha⁻¹ in páramo peatlands, making them some of the highest levels measured in mountain ecosystems (Hofstede, Segarra, and Mena 2003; Hribljan et al. 2015; Hribljan, Suarez, Heckman et al. 2016). These large concentrations of C in páramo soils explain the high levels of moisture (Farley, Kelly, and Hofstede 2004), high water retention capacity, and high hydraulic conductivity (Buytaert et al. 2006) that make them so important in the regulation and provisioning of water. Carbon storage is also clearly important in terms of climate change abatement because three times more C is stored in the soil than aboveground biomass and is a relatively stable carbon sink, especially in the deeper profile (Fontaine et al. 2007; Masciandaro et al. 2018).

Nevertheless, though soil integrity is vital for the maintenance of the health of páramo ecosystems and their ability to provide ecosystem services, it is threatened by land use change, especially by agriculture, which has been increasing across the region (Ross, Fildes, and Millington 2017). Agriculture has been shown to have negative effects on páramo soil by irreversibly changing soil properties and decreasing soil C stocks (Podwojewski et al. 2002; Farley et al. 2013), but the rate at which soil C is lost is poorly known. Traditional agricultural practices vary throughout the Andes, but farmers often follow a similar fallow-based system (Sarmiento, Monasterio, and Montilla 1993; Pestalozzi 2000). Fields are cultivated for several years, left fallow briefly, and then become pasture for livestock for nine to eleven years in order to restore soil fertility (Sarmiento, Monasterio, and Montilla 1993; Pestalozzi 2000; Sarmiento, Smith, and Monasterio 2002; Sarmiento et al. 2003). Each stage of the fallow agriculture system hosts different plant species and is subjected to very different management practices, so each stage can have varying effects on soil properties and should be considered separately when studying the effects of land use change.

Tillage during the cultivation stage, for example, has been shown to increase erosion and decrease infiltration rates in the páramo, which has negative effects on water regulation (Poulenard et al. 2001). Agriculture threatens the innate high C content, porous structure, and high hydraulic conductivity that prevent surface runoff and help direct water to subsurface catchments (Buytaert et al. 2006). Fallow periods after cultivation have a disproportionate share of water runoff and soil erosion from rainfall compared to other land uses, likely due to the lower soil organic matter (SOM) content found in these fields, which is important in moisture retention (Harden 1991, 1996). As for the effects of pastures, in puna ecosystems (dry grasslands similar to the páramo in central and southern South America) several studies have found little to no difference in soil organic carbon (SOC) stocks between pastures and natural sites (Zimmermann et al. 2010; Sylvester et al. 2017). However, Podwojewski et al. (2002) found that pasture usage in páramo ecosystems decreases water provisioning services because of losses of soil macrostructure, the development of a water-repellent microstructure, and decreases in soil C. Thus, understanding the effects of each stage of the traditional agriculture system is important in order to determine the overall effects on soil functioning because one taken by itself may offer an incomplete picture.

Despite the potential importance of agricultural practices on C storage and dynamics, there is still little information about soil C storage in the region (Scharlemann et al. 2014) and the changes that this reservoir might be experiencing under current patterns of land use change. In this context, we studied the páramo landscape surrounding the city of Quito, Ecuador, as a case study to understand the patterns of land use change in the high Northern Andes and their impact on soil C storage. Specifically, in order to assess the current state of C storage under different land use regimes, we employed a chronosequence of agricultural sites to determine rates of soil C change. This was coupled with land use change modeling to ascertain how much natural páramo land has been converted over the last thirty years and what that implies in terms of soil C storage. We hypothesized that agricultural practices would lead to C losses and that we would see a small but steady loss of C over time. Furthermore, we hypothesized that land use conversion from natural to agricultural uses has increased over the last few decades translating to large soil C losses from the Ecuadorian páramo.

Materials and methods

Study area

Soils were sampled near Quito, Ecuador, between the regions of Antisana Ecological Reserve and Cayambe–Coca National Park (Figure 1, Supplementary Materials 1). This study area is dominated by a matrix of tussock grasslands dominated by Calamagrostis intermedia and Festuca spp., interspersed with a layer of shrubs and forbs (e.g., Pentacalia, Diplostephium, Hypericum), and occasional patches of high-elevation forests (Polylepis spp. and Gynoxis spp.). Because most of our study area lies within a national park and ecological reserve, most of the páramo in and around our study area has been protected as it is the main source of water for Quito. Mean annual temperature in this area is 8.2°C with little seasonal variation, and mean precipitation recorded at 4,100 m reached 1,375 mm. Soils in this area are typical Andosols developed on volcanic ashes and with high contents of organic matter (Fehse et al., 2002).

Figure 1. Study area near Quito, Ecuador. White lines delineate areas above 3,500 m used in land use and land cover classification and subsequent C stock calculations. Black rectangles show sites of soil sampling. Rectangle S represents sampling site Sincholagua, A is sampling site Antisana, and P is sampling site Piedras Blancas. Map projection is WGS 1984/UTM 17  N.

Agriculture in the area occurs mainly between 3,500 and 4,000 m a.s.l., so our sampling was limited to this lower range of the paramo, which extends from 3,500 to 5,000 m.a.s.l. We had 3 sampling areas with 9 to 12 plots at each area, for a total of 33 plots (Figure 1, Supplementary Materials 2). A plot was defined as a single land unit uniformly and unequivocally covered by one of the land cover types defined for the study; either agriculture (i.e., crop), pasture, fallow, or natural vegetation (i.e., tussock grassland or tussock shrubland). Crop, pasture, and fallow sites are all considered agricultural uses because they are all part of the agricultural cycle practiced in this region. Natural vegetation categories represent the end points of a gradient from tussock grasslands clearly dominated by Calamagrostis intermedia to tussock shrublands in which the tussock grassland exhibits a prominent layer of shrubs. In the end, we only had two tussock grassland sites and seven tussock shrubland sites, so we decided to combine the two categories into one, hereby referred to as natural vegetation, because we did not expect significant differences between the types. For the purpose of this study, “natural” vegetation was defined as characteristic páramo vegetation that showed no sign of recent human intervention. Some of these areas, however, may have experienced burning and grazing or other management in the past. Farmers and landowners were asked about agricultural practices in the region and length of land use on the plots via informal interviews (Supplementary Materials 2). For recent human intervention estimates in the “natural,” nonmanaged sites, we relied on information from landowners who know about the history of the nearby natural sites as well as expert knowledge of the area and landscapes. We realize that in the past these sites may have been burned, grazed, or managed in some way, but by visually inspecting the plant communities in the site we determined that it had not happened recently.

Agricultural sites were all within 1 km of sites of natural vegetation used as controls to assure that there were no major differences in soil types or land structure. To control for possible differences in soil characteristics, paired sites (agricultural and natural) were selected within three distinct regions that are homogenous in their geologic origin had similar soil types throughout. The potential soil differences between regions were controlled for in later statistical analyses.

Informal discussions with landowners and workers from four different haciendas gave insight into agricultural cycles utilized in northeastern Ecuador. All haciendas maintained a crop–fallow–pasture rotation, with slight variations in the length of time of each period. Crops, almost exclusively varieties of potato (Solanum tuberosum) and fava beans (Vicia faba), are cultivated for two to three years before the area is harvested for the final time and then left fallow. The fallow field is left for a few months or up to a year. In the fallow period the field is often colonized by invasive species such as Rumex acetosella and Anthoxanthum odoratum and early successional species like Lachemilla orbiculata. After this, the field is planted with grasses, such as Lolium multiflorum and Poa annua, and continues to be enriched with species like Trifolium repens to create an improved pasture for livestock, often cows and sheep. The pasture is left for six to nine years to aid in soil fertility recovery before the cycle begins again. Owing to the almost exclusively steep slopes in the region, large machinery is not very common and much of the planting and harvesting is done with smaller machinery or by hand. Though the terrain is steep and subject to constant rain, terrace cropping is not practiced, soil erosion was evident in some sites.

Soil sampling

One deep soil core (0–5, 5–10, 10–20, 30–40, and 40–50 cm) was taken at the center of each plot with a 6.35-cm-diameter open-faced gouge auger for both carbon stocks and bulk density measurements. Páramo soils have a shallow O horizons (<5 cm), but this varies widely between sites, so this organic layer was removed before soil sampling to allow for a more equal comparison of soils between sites. Additionally, five shallow samples (0–5, 5–10 cm) were also taken with 5-cm-radius steel rings at the center of the plot and the four corners (10 m away from the center). Soil depth was measured at three to four locations on each plot and an average soil depth was calculated to calculate total C storage for the entire soil profile. Soil depth was measured using thin steel probe extensions that were hammered into the ground until the bedrock was hit.

Carbon analysis

Soil samples were dried in an oven at 60°C for four days and weighed for dry bulk density calculations (g·cm⁻³). Dried soils were then homogenized with a mortar and pestle until passed through a 2-mm sieve. Total soil C was determined via loss on ignition (LOI), which has been previously established to be as effective as measuring total C with an elemental analyzer because these soils have no carbonates (Tonneijck et al. 2010). Additionally, Hribljan, Suarez, Bourgeau-Chavez et al. (2016) found a strong correlation (R² = 0.989, p < 0.001) between LOI C content and total C measured with an elemental analyzer, indicating a similar relationship. This equation was developed by analyzing over 400 samples from our study region for total C with an elemental analyzer as well as estimating total C with LOI. We used the equation from Hribljan, Suarez, Heckman, et al. (2016) to estimate soil C in our samples because it was specifically calibrated from soils from our study site.

Samples were ashed in muffler oven for 4 hours at 500°C to determine C lost on ignition. C stocks per area (MgC·ha⁻¹) was calculated with the following equation:

$MgC\cdot h{a}^{-1}=C\cdot b\cdot d,$

where b is the sample bulk density, d is the depth of the sample, and C is the sample %C. Total C storage to the bedrock (MgC·ha⁻¹) was calculated through a partial coring approach as suggested by Chimner et al. (2014). First, we calculated the C storage for all soil depths individually. We then used the 30 to 50 cm means to estimate storage C for the layer between 20 cm and the mean soil depth for each site, which was, on average, 2.5 m. Based on previous work in the páramo in our study area, we have found that most changes in soil C in páramo soils are primarily restricted to the shallower depths (<20 cm), with little change in the deeper depths. We decided to average the deeper layers (20–50 cm) and extrapolate to the soil depth for that reason because we had fewer data points from the deeper depths. Finally, we added these four C storage values (0–5-, 5–10-, 10–20-, and 20-cm mean soil depth) to estimate total soil carbon storage for each core. Several sites only had soil samples taken from 0 to 10 cm rather than 0 to 5 and 5 to 10 cm, so these samples were excluded from the bulk density and C stock analysis because they were not comparable.

Because soils in the region are very dense, sampling with the open-face auger could alter bulk density, thus affecting C calculations. In order to assess this possibility, we compared bulk density values obtained with the auger and with steel rings in the same plots and across different land uses and found no differences (paired t test, t = −0.558, df = 92, p = 0.578) between the two bulk density methods (ring versus soil auger).

Statistical analysis

We determined the rate of C loss due to agriculture with linear regression comparing site %C to the total time of land use. Differences in %C, C stocks, and bulk density between land use types and soil depth were analyzed using linear mixed effect models with the packages “lme4” and “lmerTest” in R (version 3.5.1, R Core Team 2020). Land use and depth were set as fixed effects and region and elevation as random effects to account for spatial heterogeneity within the study area. Further analyses of differences were measured with Tukey’s post hoc tests. All data were analyzed using R statistical software.

Land use and land cover change analysis

Landsat TM (Thematic Mapper), ETM+ (Enhanced Thematic Mapper), and OLI (Operational Land Imager) images were used as data source for land use and land cover (LULC) analysis. Images from 1991, 2002, and 2017 (path 10, row 60) were selected because they were mostly cloud free and showed decadal changes in agricultural use (Supplementary Materials 4). Previously atmospherically corrected images were downloaded from the United States Geological Survey (USGS) Earth Explorer.

Maximum likelihood supervised classification was performed on all images to determine land cover of tussock shrubland, tussock grassland, bare soil, ice, developed/urban, crops, and pasture. Fallow sites were not classified because they are very heterogeneous and could easily be misclassified as bare soil or pasture. Approximately thirty training areas, with a minimum of 50 to 100 pixels each, were created per land cover type based on ground control points and expert image analysis. Masks combining any clouds in the region as well as areas lower than 3,500 m were used to clip unwanted pixels from the analysis. Peatlands are also an important ecosystem in the páramo but are only mappable with radar data (Hribljan, Suarez, Bourgeau-Chavez et al. 2016), so we excluded them from our classification. Instead, because peatland areas expand on a very slow scale and we would expect little change within twenty-six years, previously mapped peatlands (Hribljan, Suarez, Bourgeau-Chavez et al. 2016) were merged onto our classifications in order to provide a final, more accurate ecosystem map.

Accuracy of the classifications was carried out with a confusion matrix. Approximately fifty randomly selected accuracy points per class were used to determine classification accuracy. Points were validated for accuracy based on ground truthing and expert image analysis of high-resolution images. All classification and accuracy assessments were carried out using ENVI software version 4.4 (2018).

Classifications were imported into ArcGIS v.10.4 (ESRI 2002) for analysis. Areas with clouds were occasionally misclassified as urban areas due to a similar reflectance and hence were reclassified. We computed total area per land cover use with and without the area covered by clouds. Additionally, we created a transition pathway matrix in order to determine the quantity of natural tussock grassland or shrubland that had been converted into a modified land use.

Results

Carbon analysis

Total C storage for the full soil profile varied between 239 and 2,778 MgC·ha⁻¹ (Supplementary Materials 2) with a mean of 1,418 MgC·ha⁻¹ averaged across all land uses. As expected, %C significantly decreased with soil profile depth (F = 44.81, df = 4, p < 0.0001; Figure 2). There were also significant differences in %C between land use types (F = 4.67, df = 3, p = 0.008; Figure 2). Natural plots had the highest %C (12.71 ± 0.68 percent), but only between 0 and 10 cm. Beyond this depth, the %C was similar among land use types. Cultivated plots had lower %C than all other land covers (p = 0.002; Supplementary Materials 3). There was a significant negative relationship between time since initial cultivation of a plot and %C across all depths (R² = 0.11, p < 0.001; Figure 3). This indicates an absolute loss of 0.045 percent C per year due to agricultural use or a 28 percent decrease in soil C over fifty years. There was a significant negative relationship between soil C and site use age in both the 0 to 5 and 5 to 10 cm depths (R² = 0.1, p < 0.001, y = 11.3 − 0.056x and R² = 0.23, p < 0.001, y = 9.97 − 0.049x, respectively), whereas there was no significant patterns in the lower depths: 10 to 20 cm (p = 0.26, Supplementary Materials 3), 30 to 40 cm (p = 0.80), and 40 to 50 cm (p = 0.38).

Figure 2. Average soil C content for each land use type and soil profile depth. Error bars represent ±1 SE.

Figure 3. Relationship between time since initial cultivation and soil %C for all land use types.

Bulk density significantly increased with soil depth (F = 45.14, df = 4, p < 0.001), but the extent to which it increased varied among land uses (F = 4.18, df = 3, p = 0.01; Figure 4). Overall, natural sites had lower bulk density than the other land uses, including crops (p = 0.03, Supplementary Materials 3) and fallow (p = 0.02) but not pasture sites (p = 0.21). Carbon stocks (MgC·ha⁻¹) were significantly affected by depth (F = 189.49, df = 4, p < 0.001; Figure 4); however, there were also significant differences between land uses (F = 12.04, df = 3, p < 0.001). Crop sites had lower C stocks than fallow (p = 0.001, Supplementary Materials 3), pasture (p = 0.005), and natural (p = 0.03) sites, but C stocks were equal between natural and fallow (p = 0.07) and natural and pasture (p = 0.99) sites.

Figure 4. Average bulk density and C stocks for each land use type and soil profile depth. Error bars represent ±1 SE.

Land use and land cover change analysis

The land use classification was approximately 84 percent accurate with a kappa value of 0.81 (Table 1) and covered 271,545 ha. For all years, natural páramo comprised the majority of the study area, with tussock grasslands being the primary land cover (Table 2, Figure 5). As for modified land, crops or pasture were generally more common than developed land given that urban developments are rare at such elevations. Bare soil was categorized into a natural land cover because it is mainly due to lava flows around volcano tops. Bare soil increased sharply in 2002, but this was determined to be the result of a large fire in the region rather than changes in land use (Argüello 2014).

Table 1. Confusion matrix of land use classification accuracy assessment

	Ice	Bare	Water	Crops	Tussock grassland	Tussock shrubland	Pasture	Total	User’s accuracy
Ice	46	0	0	0	0	0	0	46	1
Bare	0	45	0	0	1	0	0	46	0.98
Water	0	0	46	0	0	0	0	46	1
Crops	0	0	0	27	8	7	4	46	0.59
Tussock grassland	0	0	0	0	40	6	0	46	0.87
Tussock shrubland	0	0	0	0	3	43	0	46	0.93
Pasture	0	0	0	0	19	2	25	46	0.54
Total	46	45	46	27	71	58	29	322
Producer’s accuracy	1	1	1	1	0.56	0.74	0.86

Table 2. Percentage of study area in each LULC category in 1991, 2002, and 2017

Land use	1991	2002	2017
Pasture	3.59	3.41	3.23
Crops	0.86	1.41	8.09
Developed	1.23	1.58	2.13
All modified land	5.68	6.4	13.45
Land cover
Tussock grassland	52.68	44.62	49.32
Shrubland	13.02	14.92	11.26
Bare	9.15	13.48	5.91
Ice	1.33	1.44	0.95
Grass peatland	12.93	13.39	13.57
Sedge peatland	2.80	2.90	2.94
Cushion peatland	2.11	2.18	2.21
All natural land	94.32	93.30	86.56

Note. Percentages are based on all cloud-free pixels.

Figure 5. Land cover and land use classification of páramo surrounding Quito, Ecuador, in 1991, 2002, and 2017.

Between 1991 and 2017, agricultural areas (pastures and cultivated fields) increased from 4.45 percent of the study area in 1991 to 11.32 percent in 2017, especially around the areas closest to Quito as well as to the west of the Cotopaxi volcano. Land used for cultivation increased from approximately 0.86 percent of the total study area to 8.09 percent of the total study area around Quito. This equates to an increase in agricultural land use of 838 percent over twenty-six years. However, we found an overall 10 percent decrease in pastures. The land use transition analysis (Table 3) indicated that between 1991 and 2017, 12.58 percent of originally tussock grassland was converted to agricultural use, as was 5.8 percent of tussock shrublands, signifying a 19,978-ha loss of natural páramo landscape. The conversion of areas to pastures was similarly paced across time periods, whereas land for crops was converted slowly between 1991 and 2002 but increased quickly between 2002 and 2017.

Table 3. Transition pathways from natural vegetation types to agricultural uses

	1991–2002		2002–2017		1991–2017
Pathway	Ha	%	Ha	%	Ha	%
Tussock grassland into crops	1,832.58	1.29	10,194.39	8.75	13,518.9	9.48
Tussock grassland into pasture	3,758.13	2.63	3,465.18	2.97	4,413.96	3.10
Tussock shrubland into crops	55.44	0.16	3,282.21	8.43	2,033.64	5.77
Tussock shrubland into pasture	34.02	0.09	84.6	0.22	11.88	0.03

Note. Changes indicate area converted from the natural páramo vegetation to an agricultural use.

Discussion

Total C storage was high with an average of 561 MgC·ha⁻¹ up to 1 m depths. This is high compared to other mountainous ecosystems such as the Spanish Pyrenees (59–299 MgC·ha⁻¹), the Rocky Mountains (40–220 MgC·ha⁻¹), or the Scottish Highlands (115–498 MgC·ha⁻¹; Ward et al. 2014). Other, non-peatland, páramos in Ecuador have similar C stocks; Tonneijck et al. (2010) reported 530 MgC·ha⁻¹ in northern Ecuador. In Colombian páramos, reported C levels vary from 239 to 479 MgC·ha⁻¹ and up to 530 in Venezuelan páramos (Hofstede, Segarra, and Mena 2003), indicating that the soils in our study sites have extremely high C holding capabilities. Our study also demonstrates that Andean agriculture in northern Ecuador is leading to steady losses of soil C, thus affecting dependent ecosystem services including global climate regulation (Masciandaro et al. 2018), water regulation (Farley, Kelly, and Hofstede 2004; Buytaert et al. 2006), nutrient regulation (Milne et al. 2015), food provisioning (Lal 2004), and erosion regulation (Pinheiro, Pereira, and Anjos 2004).

Carbon losses over time

By utilizing a chronosequence study with sites of different ages, we determined that a transition from the natural páramo vegetation to the crop–fallow–pasture agricultural system led to an average absolute decrease of 0.045 percent in soil C per year. Though there are no other studies in the páramo estimating rates of C loss, others have measured significant decreases in soil C from Andean agriculture. Farley et al. (2013) found the lowest soil C levels in sites previously used for agriculture rather than pastures or pine plantations. Podwojewski et al. (2002) found that overgrazing severely decreased soil C stocks of natural páramo. However, findings from studies looking at soil carbon responses to traditional Andean agriculture have been mixed. Both Abreu, Llambí, and Sarmiento (2009) and Henry et al. (2013) did not find noticeable carbon losses due to cultivation, but it was suspected that significant C was lost due to erosion. Such varied findings could be due to the fact that páramo soils are very heterogenous, and even significant losses of soil C might be masked by the extremely large pool of carbon stored in deep páramo soils and could be difficult to detect without extensive sampling (Abreu, Llambí, and Sarmiento 2009). Nevertheless, soil heterogeneity is difficult to overcome even with extensive sampling, as shown by the low, yet significant, R² value of the relationship between time since initial cultivation and soil C found in this study.

We found little difference in soil C between pasture and natural páramo sites. Natural sites had a higher %C content in the 0 to 10 cm depth, but in terms of C stocks, natural and pasture sites were similar. This was attributed to the higher bulk density in pastures due to soil compaction from livestock. Pastures have also been found to be relatively good at retaining soil C in comparison to potato fields because the extensive plant cover in pastures prevents soil erosion compared to cultivated fields (Otero et al. 2011). One study found that pastures may have SOC stocks similar to those of natural sites because of higher levels of soil weathering in pastures, which helps form organo-mineral complexes that stabilize high SOM levels even though pastures have lower plant cover than natural páramos (Heitkamp et al. 2014). Nevertheless, though C stocks may be similar between natural sites and pastures, other ESS, like water regulation, may still be negatively impacted by changes in soil structure and soil water retention found after extensive grazing (Podwojewski et al. 2002).

Fallow sites were similar to cultivated sites in terms of bulk density but had higher C levels indicating soil restoration during these agricultural rest periods. Fallow–pasture periods in highland Bolivia have been shown to completely recover soil fertility (in terms of C, N, and P) within nine to twenty-one years of fallow age (Pestalozzi 2000; Sarmiento and Bottner 2002) and others have found little to no difference in soil C content between grazed and natural puna grasslands (Gibbon et al. 2010; Zimmermann et al. 2010; Heitkamp et al. 2014). This indicates that longer fallow periods as practiced in traditional Andean agriculture may help counter some C loss caused by cultivation and should be retained as a good conservation practice. In comparison, fallow–pasture periods at our sites averaged six to nine years in length, several years shorter than the traditional cycle length. This more intensive system is likely driving the long-term decrease in soil C. Nevertheless, soil C loss seems to be primarily constrained to the upper layers of soil (0–20 cm) and the carbon deeper in the profile appears to be well protected, which may provide resilience to the páramo’s ability to store soil C.

Land use and land cover change

With our LULC analysis of the region, we were able to determine the extent of agricultural expansion around the Quito metropolitan area and its effects on C storage over the last twenty-six years. The classification had an overall 84 percent accuracy, with tussock shrublands and grasslands being the land cover classes predicted less accurately by the classification, in terms of producer’s accuracy. Though dominated by shrubs and scattered trees, tussock shrublands are also covered in tussock grasses, indicating that the reflective signature is heterogeneous between the two and will require more sampling in the future for more accurate mapping. Areas used for cultivation increased by 838 percent, or 18,511 ha, over twenty-six years, especially in the areas already close to developed zones. This expansion is not surprising because 75 percent of Ecuadorian tussock grasslands have already been converted to other land uses (Hofstede et al. 2002).

We found cropland expansion primarily in the páramo to the southeast of Quito and to the west of the Cotopaxi volcano. Growth in these areas is likely due to the proximity of already established urban areas and infrastructure as deforestation and land cover change is often associated with roads (Chomitz and Gray 1996; Nelson and Hellerstein 1997). A large highway, part of the Pan-American highway, already exists in this area and passes through the buffer zone of Cotopaxi National Park, making encroachment in the region easier. New urban development since 2003 has moved east of Quito and almost all land designated suitable for further development in Quito belongs to the eastern side as well (United Nations Environment Programme [UNEP] 2011), indicating that agricultural expansion is likely trailing urban expansion.

Pastures, in comparison to cultivated fields, experienced a 10 percent decrease in area in twenty-six years. This was unexpected given the traditionally tightly linked cycle between crops, pastures, and fallow period. Pasture usage is not increasing at the same rate as crop land, which may reflect changing agricultural practices in the region. Crissman, Antle, and Capalbo (1998) found that intensive input use in the páramo of northern Ecuador shortened the length of time required to restore soil fertility and thus the number of years that land must be left as pasture. Ecuador’s agricultural policies may also explain the decreasing pasture land. The land reform of 1954 reallocated unused arable land to people who promised to develop or cultivate it (Ross, Fildes, and Millington 2017). The average reallocated plots size was only 4.8 ha, whereas livestock require very large tracts of land (Ross, Fildes, and Millington 2017). Thus, this reallocation was found to decrease pasture size and in general promote cultivation of crops over livestock (Gondard 1988).

Future agricultural expansion is likely because Quito is growing rapidly. Between 1950 and 2010 the population increased sevenfold and it is expected to continue increasing (UNEP 2011). Agriculture is also an important job sector for the country (Egáz Yerovi 2014), signifying a possible increase in agriculture-related jobs. We found a much higher rate of agriculture expansion between 2002 and 2017 than between 1991 and 2002, suggesting that high rates of encroachment could continue into the future as well. This effect of a growing population on the páramo may also be exacerbated by climate change because a warming climate means that alpine ecosystems are moving upwards (Pauli, Gottfried, and Grabherr 1996, 2003; Kelly and Goulden 2008). It was predicted that temperatures in the páramo will rise and precipitation levels will fall (Ruiz, Martinson, and Vergara 2012), making land that was previously inhospitable now available for development. In conjunction, climate change also threatens C storage because elevated temperatures will lead to a stimulation of soil microbes (Heimann and Reichstein 2008), leading to a net decrease in SOC (Buytaert, Cuesta-Camacho, and Tobón 2011; UNEP 2011), which will further exacerbate any climate change effects.

Though we found clear signals of C loss in our study sites, likely due to lower belowground inputs from crops compared to natural vegetation and exposed soil that can lead to warming and SOM decomposition, these results may only provide conservative estimates of C losses because we did not measure C losses due to soil erosion. Soil erosion was very evident at many of our sites, because the highly sloped sites facilitate soil loss when plant cover is low. Whether or not soil erosion is an ultimate C sink or source is highly contested (Jacinthe and Lal 2001; Lal 2005; Van Oost et al. 2007), but we argue that in our system, soil erosion is likely a C source. Crop cultivation leads to large amounts of bare soil in our sites, which, upon exposure to the elements, can lead to soil drying, which irreversibly affects Andosol structure (Buytaert et al. 2002). Moreover, there are several heavy rainfall events each season that can wash away bare soil in cultivated or fallow fields on steep slopes, washing the soil to lower elevations. Thus, it is possible that C losses from erosion may be a significant source of C emissions in the páramo, but there is little to no research currently available. Future studies looking at the rates of C loss due to agricultural erosion in the region, from either wind or rainfall, would help provide a more accurate estimate of the changing C budget in the region.

Conclusion

The cultivation of crops has detrimental effects on C stocks in the Ecuadorian páramo. This study has provided the first estimate of C loss due to agricultural practices in the highlands of Ecuador and can be a starting point for more effective páramo conservation efforts in the region by determining the regions at most risk for further land use change. Land use change has increased rapidly in the last decade and is also likely to be accelerated by climate change, indicating even more potential C loss, which will negatively affect both regional and global ecosystem services. However, because agriculture is such an important livelihood for a large part of the Ecuadorian population, local management must consider ecosystem-based adaptation measures to help support communities. Further studies looking at the efficacy of agricultural methods that reduce soil and C loss, such as traditional longer fallow and pasture periods, should be explored as possible alternatives to prevent further soil C loss.

Acknowledgments

Publication made possible in part by support from the Berkeley Research Impact Initiative (BRII) sponsored by the UC Berkeley Library. We thank the Freshwater Ecology lab at USFQ for the field and laboratory space and equipment. We also thank Shahar Amitay and Valerie Dong for their help in the field and laboratory and Timothy Bowles for helpful comments. We also extend our gratitude to the U.S. Forest Service and Sustainable Wetlands Adaptation and Mitigation Program for use of the peatland data and maps. Finally, we are thankful for the anonymous reviewers for their comments on the article.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental material for this article can be accessed on the publisher’s website.

Additional information

Funding

This study was financially supported by the European Commission through the program Erasmus Mundus Master Course—International Master in Applied Ecology (EMMC‐IMAE; FPA2023‐0224/532524-1-FR-2012-1-ERA MUNDUS-EMMC). Support for ES and SC was provided by a Research Grant provided by Universidad San Francisco de Quito.