Forest cover change and its carbon dynamic of the karst area in Bulusaraung, South Sulawesi, Indonesia

Abstract

Information on land forest cover changes and carbon stocks in karst areas is vital for management planning. Therefore, an analysis of the dynamics of land forest cover, carbon stocks of secondary karst forests, and their future projections is necessary. This study aimed to analyze the carbon stocks of secondary forests, forest cover change, and their projections in karst areas. We used 43 sample plots of 20 m x 20 m to measure Above Ground Biomass (AGB) and soil organic carbon stocks. Forest cover data were collected using remote sensing and GIS tools. Data analysis of the ABG of the secondary forest for the karst plain and karst tower was conducted using the Chave equation. Land-cover change was analyzed by delineating the 2011, 2016, and 2021 land-cover data from satellite imagery and testing the accuracy using survey data fields and high-resolution images. Land-cover projections for 2026 and 2031 were obtained using Cellular Automata Markov (CA-Markov) analysis. The results revealed ten land cover types in the karst area of the Bulusaraung Forest Management Unit (FMU). Secondary dryland karst forests (tower and plain karst) were the dominant forest types (60%) in the Maros-Pangkep essential karst ecosystem area. Both karst types of secondary forests are relatively high carbon stores. Data projection shows that the mining area will grow extensively by 2031. Secondary dryland karst forests will be degraded by 15.5% in 20 years due to increased mining activities, conversion to paddy fields, dryland agriculture, and mixed dryland agriculture. Karst ecosystems are vital for carbon and water storage. Therefore, some strategies are required for the sustainable management of Maros-Pangkep Karst are as follows: (i) controlling the opening of mining areas according to sustainable environmental principles, (ii) encouraging the implementation of post-mining land rehabilitation, (iii) increasing public awareness of the importance of preserving karst areas, and (iv) implementing soil conservation rules by prioritizing agroforestry systems in dryland farming.

Keywords:

• carbon sequestration
• secondary karst forest carbon
• forest cover change
• Maros-Pangkep Karst
• UNESCO Global Geopark

1. Introduction

In recent decades, karst areas have received special attention from experts for their ability to mitigate climate change by storing carbon (C) pools and reducing atmospheric CO₂. Karst areas account for 15% of the earth’s landmass (Song et al., 2017), and approximately 25% of the world’s population depends on ecosystem services (Xiao & Weng, 2007). Karst areas are important carbon sinks in the terrestrial pool and play a significant role in the global carbon cycle and climate change. Zhang et al. (2022) stated that karst areas have a significant carbon sink capacity and contribute 29.4% of the terrestrial CO₂. Falkowski et al. (2000) reported that karst areas had more than 6 × 107 Pg C. In other words, 99% of the total carbon on earth is stored in carbonate rocks (Falkowski et al., 2000), which constitute the largest carbon stores on earth (Cao et al., 2016; Martin, 2017). Carbonate rocks also store carbon for millions of years (Wu & Wu, 2022).

In Indonesia, 154,000 km² of tropical karst areas streches from Sumatra to Papua (Samodra, 2003) with high rainfall, and intensive karstification processes significantly contribute to carbon storage. Haryono (2011) stated that in Indonesian tropical karsts, every 1 ton dissolution of limestone is followed by the absorption of 0.12 tons of CO₂ from the atmosphere. For example, Jonggrangan karst absorbs 12,900 tons of CO₂ yr⁻¹, while Gunungsewu karst absorbs 293,800 tons of CO₂ yr⁻¹ (Cahyadi, 2010).

The carbon storage process in karst areas involves not only non-living elements and hydrogeological processes, but also biotic elements and biological processes. Cao et al. (2023); Wu and Wu (2022); Yang et al. (2015) reported that photosynthetic organisms play vital roles in carbon sink processes. Karst forest ecosystems are widely recognized as crucial agents for absorbing greenhouse gases and mitigating climate change (Kang et al., 2020; Kang & He, 2013). Various vegetation types that form karst forests play essential roles in carbon sequestration (Huang et al., 2015; Liu et al., 2020; Ni et al., 2015; Wu & Wu, 2022) and enhance the carbon sink capacity of ecosystems (Yang et al., 2022). The capability and stability of karst carbon sinks increase significantly with the presence of forest cover (Jackson et al., 2005; Tong et al., 2018). Li (2012) and Luo et al. (2005) stated that karst vegetation restoration increased karst carbon sinks by 19–23%. Furthermore, Zhang (2011) stated that land-cover transformation in karst areas affects carbon sinks. Karst carbon sinks in secondary and primary forests were three and nine times higher, respectively, than those in shrubland.

However, research on the roles of karst forests in carbon storage in Indonesia is limited (Danardono et al., 2019; Sari et al., 2021). A lack of comprehensive research has led to inadequate data on carbon storage in national-level karst forests. This information is vital for climate change mitigation, biodiversity conservation, fulfilling international reporting obligations, and guiding evidence-based policymaking and resource management decisions. An understanding of carbon storage not only facilitates progress tracking and impact assessment of land-use changes but also contributes significantly to global commitment in combating climate change. Therefore, this study analyzed the role of forests in the Maros-Pangkep karst in carbon storage and sequestration. This study was conducted in the Maros-Pangkep karst area, which is essential for Indonesia and the world. The Maros-Pangkep Karst is part of the Maros-Pangkep UNESCO Global Geopark, one of 195 UNESCO global parks in 48 countries (International Geoscience and Geoparks Programme, 2023). This tower karst is higher than 200 m (Srijono, 2008) and is the only tower karst in Indonesia (Yusran, 2017).

Research on the roles of forests in the Maros-Pangkep karst area in carbon storage is essential because Drew (2017) and Zhao et al. (2020) state that karst is a fragile ecosystem and is easily damaged by human activities. Additionally, the pressure and negative impacts caused by human activity in karst areas are relatively high (Jiang et al., 2014; Parise & Gunn, 2007) and can disrupt and change karst areas. This condition was also observed in the Maros-Pangkep karst area. In this area, the potential for environmental change, forest degradation, and changes in carbon storage due to human activity is high. This can be seen from the activities of local communities that change forests for agriculture to meet their daily needs, and even use this area for large-scale limestone mining activities, which produce cement and marble for national demand (Budiyanto et al.; Fatinaware et al., 2019). It is crucial to determine the dynamics of changes in forest cover in the essential ecosystem of the Maros-Pangkep karst. This change can be observed in vegetation cover dynamics (Peng et al., 2021). Therefore, this study aimed to (1) determine carbon storage in secondary forest karst areas in the Maros-Pangkep Karst and (2) analyze forest carbon change and its projections based on forest cover dynamics in the karst area managed by the Bulusaraung Forest Management Unit (FMU). We hypothesized that the aboveground carbon of forest is higher than other land covers, and the land cover will gradually change into more open area. Remote sensing and Geographic Information Systems (GIS) were used to identify forest cover changes and predict future land cover projections. These technologies have been recognized as effective over a certain period (Obiahu et al., 2021; Setya & Wiryani, 2019). Our findings provide recommendations for FMU managers for better and more sustainable karst ecosystem management.

2. Materials and Methods

2.1 Study area

The Maros–Pangkep karst area is a part of the Maros–Pangkep UNESCO Global Geopark. This karst area has an area of approximately 46,200 ha and covers five regencies: Pangkajene, Kepulauan (Pangkep), Maros, Barru, Bone, and Jeneponto (Achmad & Hamzah, 2016). This research was conducted in November 2021 to Januari 2022 in the Maros-Pangkep karst area, which is managed by the Bulusaraung FMU under the coordination of the South Sulawesi provincial government. The Bulusaraung FMU was established based on South Sulawesi Governor Regulation Number 134 on 2017 September 8, 2017. Bulusaraung FMU manages protected and production forests, which are administratively located in two regencies, Maros and Pangkep, South Sulawesi Province, Indonesia. The study area and karst distribution are shown in Figure 1.

Figure 1. Study area and the karst distribution

2.2 Forest carbon measurement

Based on the karst formation in the forest area, the study area was divided into two types: 1) plain karst, which is located in the Maros Regency and is characterized by small and low features spread across the study site. Most are underground, so only a few are visible; 2) tower karst, which is located in the Pangkep Regency with high and steep cliffs reaching up to tens of meters. In the tower karst, data collection was absent when the elevation was nearly 90°. We recorded data only from secondary forests located in accessible areas. A total 43 sample plots were established: 20 on the karst plain and 23 in the karst tower. The plain and tower karsts in the Maros-Pangkep area are shown in Appendix 1.

The methods for data collection and analysis were based on SNI 7724:2011 (Badan Standardisasi Nasional, 2019) for ground-based forest carbon accounting. Field data were collected from rectangular and square plots. The selection of plots depended on the field conditions. For example, rectangular plots were generated when the observation site was located at the edge of two karst cliffs. Rectangular plots were selected when they were located in narrow areas (Putri et al., 2017). Nevertheless, the minimum size of the plots remained uniform at 400 m² (Badan Standardisasi Nasional, 2019).

Vegetation data were sampled by applying a line transect method (Fachrul, 2012) at a distance of 50 m. Inside of each plot, we established 4 (four) nested sub-plots for the trees (20 m × 20 m), poles (10 m × 10 m), saplings (5 m × 5 m), and understories and seedlings (2 m × 2 m) in each plot. We defined (i) tree as woody vegetation with a diameter at the breast height (DBH) > 20 cm, (ii) pole as woody vegetation with DBH > 10 cm and < 20 cm, (iii) sapling as woody vegetation with DBH 2 cm – 10 cm, and (iv) seedling as woody vegetation with DBH < 2 cm and height ≤ 1,5 m; and (v) understory as vegetation grown on forest floor including herbs, ferns, and shrubs (Badan Standardisasi Nasional, 2019). During fieldwork, we recorded all data including the name of the species, DBH, and the total height of the plant.

We also collected data on other carbon pools, including (i) herbs or understory, (ii) litter, (iii) standing dead wood, (iv) fallen dead wood, and (v) soil organic matter. Herbs or understory data were obtained from plots of 2 m × 2 m. The aboveground biomass of the herbs and understory were cut and collected destructively. The total biomass was then weighed and 300 g was sub-sampled and stored in a plastic bag for laboratory analysis (Badan Standardisasi Nasional, 2019).

Necromass measurements were performed in 2 m × 2 m plots. All the litter samples were collected and weighed. A total of 300 g of litter was subsampled and stored in plastic bags for laboratory analysis. A geometric method was applied to measure necromass of dead standing trees within plots of 20 m × 20 m. Data on DBH and tree height were recorded in a notebook. The wood samples were collected for laboratory analysis. For fallen dead wood, necromass was measured based on the volume. We collected data in plots of 20 m × 20 m, including ground and top diameters and wood length (Badan Standardisasi Nasional, 2019).

Soil organic matter samples were taken from 0 cm – 10 cm depth in plot of 20 m x 20 m. Samples were derived from the center and each of the four corners of the plots used destructive sampling method. All collected samples were then combined (Badan Standardisasi Nasional, 2019). C organic stock was analyzed using Walkley and Black method at the Soil Laboratory of the Agricultural Faculty of Hasanuddin University. We also used soil ring samples to collect soil from 0 cm – 10 cm depth to measure the bulk density of soil.

Aboveground biomass in the secondary forest of karst plain and tower karst was determined using an allometric equation (1)(1) $$\mathit{\text{AGB}}=0.0509\hspace{0.17em}x\hspace{0.17em}\rho \hspace{0.17em}x\hspace{0.17em}{D}^2\hspace{0.17em}x\hspace{0.17em}H$$(1) developed by Chave et al. (2005). (1) $$\mathit{\text{AGB}}=0.0509\hspace{0.17em}x\hspace{0.17em}\rho \hspace{0.17em}x\hspace{0.17em}{D}^2\hspace{0.17em}x\hspace{0.17em}H$$(1) where:

AGB = aboveground biomass (kg tree⁻¹)

D = diameter at breast height (cm)

ρ = wood density (g cm⁻³), derived from ICRAF database (http://db.worldagroforestry.org/wd)

H = tree height (m)

Carbon from each carbon pool, including understory, seedling, sapling, pole, tree, litter, and necromass, was analyzed using the formula: (2) $${\mathbf{C}}_{\mathbf{b}}=\mathbf{B x \%C }\mathbf{\text{organic}}$$(2) where:

Cb = Carbon of biomass (kg)

B = Total biomass (kg)

% C organic = the carbon percentage (0,47)

Carbon Dioxide (CO₂) absorption was estimated using formula (2) from Brown (1997) to calculate CO₂ absorption: (2) $${\mathbf{\text{CO}}}_2=\frac{{\mathbf{\text{MrCO}}}_2}{\mathbf{\text{Ar}}\mathrm{ }\mathbf{C}}\text{x}\mathbf{C}\mathit{b}$$(2)

Where:

CO₂ = carbon dioxide absorption

Mr CO₂ = CO₂ relative molecular weight (44)

Ar C = C molecular weight (12)

Cb = carbon biomass

An independent sample T-test was used to test the hypothesis that there are differences in the carbon stock and carbon sequestration of karst plains and karst towers.

2.3 Forest cover changes

The first stage was to determine the area of interest using a boundary map of the Bulusaraung FMU and a map of the Maros-Pangkep karst area. ArcGIS software was used to intersect the map of the Bulusaraung FMU area with the map of the Maros-Pangkep karst essential ecosystem area. The resulting map indicated the karst areas within the FMU working area. Next, Landsat image preparation and land cover delineation were performed using ArcGIS software. Landsat 7 ETM + path 114 row 63 acquisition September 21, 2011, Landsat 8 OLI/TIRS path 114 row 63 acquisition September 10, 2016, and Landsat 8 OLI/TIRS path 114 row 63 acquisition August 23, 2021, were used for land cover delineation. A visual method (onscreen digitization) was applied for image interpretation (Jensen, 2015) for the 2011, 2016, and 2021 time series. The year 2011 was chosen as the starting point for the observations considering the availability of high-resolution satellite images for accuracy testing.

The pan-sharpened mode was used with a spatial resolution of 15 m for image interpretation (Gilbertson et al., 2017). A land-cover survey was conducted to obtain a training area to determine the type of land cover. Information on land cover type was obtained directly from visual observations in the field assisted by a drone. Field data locations were recorded using a handheld global positioning system (GPS).

The land cover maps generated from the interpretation of Landsat imagery must be tested for accuracy based on the results of the field check and by verifying the appearance of high-resolution images from Google Earth Pro from September 2011 to October 2021. An accuracy test was performed using the confusion matrix method (Lillesand et al., 2015), which compared the types of cover interpretation results with the actual conditions in the form of a sample plot. The accuracy is expressed by the Kappa coefficient (K), which is calculated using Equation (3)(3) $$\mathrm{K}=\frac{\mathrm{N}\sum _{\mathrm{i}=1}^{\mathrm{I}}{\mathrm{X}}_{\mathrm{\text{ii}}}-\sum _{\mathrm{i}=1}^{\mathrm{I}}{\mathrm{X}}_{\mathrm{i}+}{\mathrm{X}}_{+\mathrm{i}}}{{\mathrm{N}}^2-\sum _{\mathrm{i}-1}^{\mathrm{I}}\left({\mathrm{X}}_{\mathrm{i}+}{\mathrm{X}}_{+\mathrm{i}}\right)}$$(3) (Lillesand et al., 2015; Richards & Jia, 2006). (3) $$\mathrm{K}=\frac{\mathrm{N}\sum _{\mathrm{i}=1}^{\mathrm{I}}{\mathrm{X}}_{\mathrm{\text{ii}}}-\sum _{\mathrm{i}=1}^{\mathrm{I}}{\mathrm{X}}_{\mathrm{i}+}{\mathrm{X}}_{+\mathrm{i}}}{{\mathrm{N}}^2-\sum _{\mathrm{i}-1}^{\mathrm{I}}\left({\mathrm{X}}_{\mathrm{i}+}{\mathrm{X}}_{+\mathrm{i}}\right)}$$(3)

Remark:

N = The total number of samples in the confusion matrix

Xii = number of samples on the main diagonal of the confusion matrix

Xi+ = the number of samples of all columns in a row

X + i = number of samples of all rows in a column

After obtaining land cover maps for 2011, 2016, and 2021, projections for 2026 and 2031 were obtained using the cellular automata (CA) system approach. A CA is a discrete dynamic system that works spatially and temporally according to a specific pattern in the same area (Hand, 2005). CA has been successfully integrated into a model in IDRISI software, making it easier to use in LULC modelling (Hadi et al., 2014).

Land cover projection was performed using the Cellular Automata Markov (CA-Markov) module in the IDRISI-Selva software (Mannan et al., 2019). Initially, the Markovian transition estimator was run with land cover raster maps as inputs for 2011 and 2016 to produce the following: 1) a Markov transition area file (a text file that contains the number of pixels that are expected to change from one land cover type), and 2) image files of transition suitability. Furthermore, the projection with iterations was set to five to produce a map of land cover projections for 2021, and the map was validated using the existing (2021) land cover map. Validation was performed in two ways: utilizing the validation and Crosstab modules on IDRISI-Selva. If the result is higher than 0.8, the land-cover projection process for 2026 and 2031 can continue. Changes in land cover and projection maps were analyzed to determine the trends of spatial and temporal changes (Ait El Haj et al., 2023; Aneesha Satya et al., 2020).

3. Results

3.1 Forest Biomass, Carbon stock, and Sequestration

The results showed differences in the volume of carbon stocks stored in vegetation, necromass (litter, standing dead wood, and fallen dead wood), and soil (Table 1). Carbon measurements of the vegetation showed that the tree stage contributed the most to biomass and carbon stocks. Carbon stock measurements conducted on the necromass showed that the litter shared the largest proportion of carbon stocks. The results also showed that the soil carbon stock in karst areas was higher than that in vegetation and necromass.

Table 1. Biomass and carbon stock of karst plain and karst tower

No	Level	Karst plain			Karst tower
		Biomass		Carbon	Biomass		Carbon
		ton/ha	%	ton/ha	ton/ha	%	ton/ha
A	Vegetation
1	Herb	0.92	0.55	0.43	0.71	0.44	0.33
2	Sapling	13.13	7.83	6.17	12.09	7.49	5.68
3	Pole	23.44	13.97	11.02	22.58	13.98	10.61
4	Tree	128.70	76.72	60.49	125.27	77.57	58.88
	Total 1	166.20	99.07	78.11	160.65	99.47	75.51
B	Necromass
5	Litter	1.55	0.92	0.730	0.8500	0.52	0.4000
6	Standing dead wood	0.011	0.01	0.005	0.0009	0.00	0.0004
7	Fallen dead wood	0.002	0.00	0.001	0.0013	0.00	0.0006
	Total 2	1.56	0.93	0.735	0.848	0.53	0.399
8	Soil			870.20			1,040.39
	Grand Total	167.76	100	949.05	161.50	100	1,116.29
	Median			76.80			62.86
	Average			78.85ns			75.90ns
	Std deviation			38.17			38.53

The results of independent sample T test analyses showed that there were no significant differences in carbon stock and carbon sequestration between the karst plains and karst towers (Table 1 and Table 2). T-test analysis is shown in the Appendix 2.

Table 2. Carbon sequestration of each type of area.

No	Type of area	Carbon sequestration (ton CO2)
		Herb	Sapling	Pole	Tree	Total
1	Karst plain Maros	1.58	22.64	40.39	221.80	286.41
2	Karst tower Pangkep	1.22	20.84	38.91	215.89	276.85

Note: all values between karst plain and karst tower are not significantly different according to T-test analysis

3.2 Forest cover changes

Land-cover delineation based on the interpretation of satellite imagery resulted in land-cover maps for 2011, 2016, and 2021 (Figure 2).

Figure 2. Map of land cover changes 2011 – 2021

Field observations of 81 points were made to obtain an overview of the land cover and were used to calculate the accuracy of the resulting map. The accuracy test result in Table 3 showed a Kappa index of 94.2% and overall accuracy of 96.3%; indicating the map resulted from the interpretation was acceptable.

Table 3. The result of accuracy test

Land cover		Field Check
		SDR	St	MA	DA	MDA	GL	PF	Sh	FP	OW	Total	User accuracy
Interpreation	SDR	43	–	–	–	–	–	–	–	–	–	43	100
	St	–	1	–	–	–	–	–	–	–	–	1	100
	MA	–	–	1	–	–	–	–	–	–	–	1	100
	DA	–	–	–	2	–	–	–	–	–	–	2	100
	MDA	1	–	–	–	11	–	–	1	–	–	13	84.6
	GL	–	–	–	–	–	1	–	–	–	–	1	100
	PF	–	–	–	–	–	–	1	–	–	–	1	100
	Sh	1	–	–	–	–	–	–	16	–	–	17	94.1
	FP	–	–	–	–	–	–	–	–	1	–	1	100
	OW	–	–	–	–	–	–	–	–	–	1	1	100
	Total	45	1	1	2	11	1	1	17	1	1	81
	Producer accuracy	95.6	100	100	100	100	100	100	94.1	100	100
	Kappa index	94.2
	Overall accuracy	96.3

Remarks: SDR = Secondary dryland forest; St = Settlement; MA = Mining area; DA = Dryland agriculture; MDA = Mixed dryland agriculture; GL = Grassland; PF = Paddy field; Sh = Shrub; FP = Fishpond; OW = Open water

Furthermore, predictions of land cover in 2021 were made based on land-cover maps for 2011 and 2016. The prediction result was validated using the IDRISI-Selva software by comparing it with the actual land cover map for 2021. The results of the validation module as shown in Fig. 3 indicates that the K-standard value is 0.93. Using the crosstab module, the validation also resulted in a very high value, as indicated by the overall Kappa Index of Agreement (KIA) of 0.93.

Figure 3. Validation results of land cover projection in 2021

Based on the feasibility of the projection model, the following predictions were made for the land cover in 2026 and 2031, as shown in Fig. 4.

Figure 4. Land cover map projections of 2026 and 2031

Based on the results of the land-cover mapping from 2011 to 2031, the area dynamics for each land-cover class are illustrated in Figure 5.

Figure 5. The dynamic area of each land cover type

Appendix 1. Two types of karst at the Maros – Pangkep area. a. Plain karst; b. Tower karst.

Figure 5 shows that dryland agriculture, mining areas, and paddy fields tend to increase in size; secondary dryland forests, mixed dryland agriculture, and shrubs tend to decrease, whereas others tend to stagnate. Secondary dryland forest is the land-cover type with the most significant decreasing trend, whereas mining is experiencing the fastest increasing trend. Details of the area and percentage of each land cover class are listed in Table 4.

Table 4. Actual and projection of land cover at the karst area of Bulusaraung FMU

	Actual land cover (ha)						Projection (ha)				2011-2031 changes	% ^(*)
	2011	%	2016	%	2021	%	2026	%	2031	%
Secondary dryland forest	4622.7	59.2	4613.1	59.1	4579.2	58.6	4222.2	54.1	3904.4	50.0	718.3	−15.5
Dryland agriculture	166.8	2.1	182.5	2.3	219.6	2.8	418.1	5.4	510.7	6.5	343.9	206.2
Mixed dryland agriculture	1111.8	14.2	1109.8	14.2	1091.1	14.0	1089.7	14.0	1075.9	13.8	35.9	−3.2
Settlement	0.9	<0.1	0.9	<0.1	1.6	<0.1	1.6	<0.1	1.6	<0.1	0.7	72.6
Mining area	47.2	0.6	56.8	0.7	93.8	1.2	277.3	3.6	506.3	6.5	459.0	971.6
Shrub	1805.5	23.1	1788.5	22.9	1751.6	22.4	1629.2	20.9	1534.9	19.7	270.5	−15.0
Grassland	5.9	0.1	6.9	0.1	7.6	0.1	7.4	0.1	7.4	0.1	1.5	24.9
Paddy field	46.4	0.6	48.7	0.6	62.7	0.8	161.9	2.1	266.1	3.4	219.7	473.4
Open water	1.6	<0.1	1.6	<0.1	1.6	<0.1	1.5	<0.1	1.5	<0.1	0.1	−4.8
Fishpond	1.0	<0.1	1.0	<0.1	1.0	<0.1	1.0	<0.1	1.0	<0.1	0.0	4.2
Total	7809.8		7809.8		7809.8		7809.8		7809.8

⁽*⁾a negative number indicates a decrease in area

Table 4 indicates a projected decrease in the area of forest cover in 2031 of 718.3 ha or a decrease of 15.5% compared to 2011. Based on the results of the spatial analysis, changes in forest cover were due to change in other land-cover types, as presented in Table 5.

Table 5. Details of changes in forest cover to/from other land cover

Land cover sources	Decrement (ha)	Increment (ha)	Change (ha)	Change proportion (%)
Dryland agriculture	104.25	0.06	104.19	14.51
Mixed dryland agriculture	83.57	0.91	82.66	11.51
Settlement	0.08	0.02	0.06	< 0.01
Mining area	316.01	0.00	316.01	44.00
Shrub	17.33	7.72	9.61	1.34
Grassland	0.16	0.42	−0.26	< 0.01
Paddy field	206.04	0.00	206.04	28.69
Open water	0.01	0.03	−0.01	< 0.01
Fishpond	0.03	0.05	−0.02	< 0.01
Total	727.48	9.20	718.27	100

Table 5 shows that changes in the area of forest cover are more dominant due to area conversion of mining, paddy fields and dryland agriculture. This condition was not offset by the increase in forest cover area, which in total only increased by 9.2 ha in a period of 20 years, mainly in the form of shrubs.

In line with Table 4, based on the predictions of the dynamics of the secondary forest area on karst and the dynamics of karst forest carbon stock, a decrease in both the total area of secondary dryland forest and carbon stock in the two types of karst topography (karst tower and karst plain) is shown in Table 6.

Table 6. Dynamic carbon stock prediction based on land cover change prediction

Years	Area of secondary dryland forest (ha)		Carbon stock (ton)
	Karst tower	Karst plain	Karst tower	Karst plain	Total
2011	3075.1	1547.6	5308303.5	2617236.6	7925540.1
2016	3067.6	1545.5	5295254.5	2613785.5	7909040.0
2021	3044.8	1534.4	5255939.7	2594971.5	7850911.2
2026	2805.7	1416.5	4843156.9	2395626.3	7238783.3
2031	2582.4	1322.0	4457771.0	2235734.9	6693505.9

4. Discussion

4.1. Karst forest carbon stock

The secondary forest on the Maros-Pangkep karst is quite good as a carbon storage. This can be seen from the slightly higher carbon stocks in the secondary karst forest in Bulusaraung FMU than in the karst secondary forest in Southwestern China, which was 70.3 Mg ha⁻¹ (Liu et al., 2016). However, the carbon stocks of the secondary forest of karst plain Maros and the secondary forest of tower karst Pangkep which are 78.11 ton ha⁻¹ and 75.51 ton ha⁻¹ respectively, are lower than that of secondary karst forest at Cirebon, which was 87.18 ton ha⁻¹ (Yulistyarini et al., 2016). Zhang et al. (2022) stated that there are differences in the amount of soil organic carbon in karst forests caused by differences in the successional stage and soil carbon input. Other factors that might affect the variation in carbon stock at different locations include soil thickness, the composition of plant species at the two locations, and the difference in the litter amount, which is part of the primary source of aboveground carbon input.

There were slight differences in the biomass, carbon stock, and carbon sequestration of secondary forest vegetation between the karst tower and the karst plain. The karst plain biomass, carbon stock, and carbon sequestration were slightly higher than those of the karst tower because of physical environmental conditions, such as topography and soil thickness. Karst rocks sporadically appear as irregular cones or small hills on the karst plain. Most karst bodies are located beneath thick layers of soil. The soil resulting from the sedimentation process over thousands of years had relatively high organic matter and nutrient content. This allows plants that live above the karst plain area to have a better growth rate than those in the karst tower; this results in the karst plain having a higher biomass, carbon stock, and carbon sequestration than the karst tower.

Soils in the secondary karst forest area of Maros play a vital role in carbon storage. The soil carbon stock content in the Maros-Pangkep karst secondary forest area (Table 1) was also higher than that carbon in several karst locations in Southwestern China, as reported by Wang et al. (2021). The karst plain contained a higher aboveground biomass and C-stock (from vegetation and necromass) than the karst tower; however, it contained a lower soil C-stock (Table 1). The carbon content variation in the soil SOC in tower karst and plain karst may be due to differences in the slope gradient. This is consistent with the findings of Zhang et al. (2019), who stated that differences in soil thickness, slope gradient, and rock exposure caused differences in SOC. In addition, community activity factors play an important role in causing differences in SOC (Ran et al., 2023). Zhang et al. (2020) stated that human disturbance is a strong driving force influencing the SOC. Li et al. (2020) stated that anthropogenic factors affect the amount of leaf litter that becomes soil carbon input. Areas with minimal human activity had more leaf litter than areas with high human activity. The sloping topography of the karst plain makes this area more accessible to the community than the karst towers.

4.2 Land cover change of the FMU Bulusaraung

Social-ecological systems are greatly influenced by interactions between humans and nature which have long-term impacts on environment sustainability. Human driving factors have significantly changed natural conditions including karst ecosystem in most regions of the world. High-intensity human activities have placed greater pressure on the natural environment to regulate itself and restore its condition. Consequently, human actions and interactions emerge as the main forces driving change and reshaping of national territorial patterns. These phenomena also accurs in the research site where human activities have significantly change the shape of karst ecosystem (Li et al., 2024). In the research site, land cover analysis for the period 2011–2021 shows ten types of land cover in the karst area of the Bulusaraung FMU. The most dominant existing land cover type is secondary dryland forest with a composition of almost 60%, followed by shrubs and mixed dryland agriculture. The other seven land cover types (Table 4) have a small area, with each accounting for less than 3%. The three smallest land cover types were settlements, open water, and fishponds. The contribution of the three was less than 0.1% in the period 2011–2031 with constant growth, and then tended to decrease.

The land cover change analysis also showed that the highest percentage of land cover reduction occurred in secondary dryland forest cover (15.5% in 20 years). Spatial analysis of land cover projections showed a decrease from 4,623 ha of secondary dryland forests in 2011 to 3,904 ha by 2031 (Table 4). Table 5 shows that 718.3 ha forest cover change was mainly caused by deforestation due to mining activities (44%), followed by conversion to paddy field (28.7%), dryland agriculture (14.5%), and mixed dryland agriculture (11.5%). This findings are in accordance with the study from Lezhava et al. (2020) in Zemo Imereti Plateau, located in the easternmost part of the limestone region of western Georgia and is Caucasus’ only karst plateau who revealed that the karst landscape has been completely transformed due to human strong economic activities, mainly related to the manganese unplanned, predatory extraction. The occurrence of deforestation and degradation in karst areas affects the role of karst forests in carbon sequestration. The decreasing trend of the secondary dryland forest cover is also occurring nationally. Government data shows the area of secondary dryland forest cover decreased from 38,979,000 ha; 36,506,100 ha; 35,960,600 ha in 2011, 2016, and 2021 respectively. However, the trend of increasing area occurred in mining areas and paddy fields nationally. Mining area increased from 539,500 ha; 677,000 ha; and 981,000 ha in 2011, 2016, and 2021, respectively. Paddy fields also escalated steeply from 7,118,200 ha; 7,841,200 ha; and 8,088,600 ha in 2011, 2016, and 2021, respectively (Mutiara et al., 2018; Rahayu, 2012; Ratnasari & Tosiani, 2022).

The land cover for paddy fields and dryland agriculture in the research site has increased extensively. Based on field observations, these agricultural lands developed on moderate to gentle slopes of hills or tower karsts. Although land conversion into agricultural land has a lower impact than mining activities (Gbedzi et al., 2022), agricultural activities closer to karst areas increase the threat to karst sustainability. A closer distance between karst and community activities increases the intensity of the community accessing karst areas and utilizing various natural resources in the karst ecosystem to support their livelihoods (Breg Valjavec et al., 2018). However, karst ecosystems are sensitive to agricultural activities due to pesticide and chemical fertilizer application as well as improper soil conservation practices. These activities will increase runoff, erosion rates, water pollution, and weaken climate change mitigation efforts (Jelev et al., 2021).

Another change in the secondary karst forest was predicted to be mixed dryland agriculture (approximately 11.51%). Although forests have been converted into mixed dryland agriculture over a period of 20 years (Figure 5 and Table 4), there has been a decline in the area of mixed dryland agriculture. This is because the rate of increase in mixed dryland agricultural areas originating from forest cover is smaller than the rate of conversion for other uses, particularly dryland agriculture. Mixed dryland agriculture is an agroforestry system that is dominated by woody plants and fruit trees. In some locations, woody plants are combined with food crops and bamboo plants. The change in forest area to mixed dryland agriculture indicates forest encroachment by communities, similar to that of another FMU in Indonesia reported by Rijal et al. (2019); Yurike et al. (2021). Based on further analysis, after this change, the area of mixed dryland agriculture is predicted to decrease by 3.2% by 2031. Mixed dryland agriculture with an abundance of long-lived trees can play a role in carbon storage and climate change mitigation (Ghale et al., 2022). Moreover, our field observations showed that the communities applied the concept of agroforestry, which is characterized by the dominance of woody plants with wood and fruit production orientations. Mixed woodlots improved soil properties and might be considered as a technique for restoring degraded dryland soils while also acting as key carbon dioxide and nitrogen sinks (Syano et al., 2023). In some locations, woody plants are combined with bamboo or other food crops. Areas with many long-lived trees, such as mixed dryland agricultural areas, can absorb and store carbon for long periods of time.

Based on the projections analyzed until 2031, the Maros-Pangkep karst area managed by the Bulusaraung FMU faces a threat of significant land use change. Our predictions show that the secondary karst forest area will decrease and be replaced by other areas, such as mining areas, which will contribute to the most significant increase. In 2011, the mining area was only 47.2 ha or 0.6% of the total area, and it is projected that the mining area will increase by up to 459 ha or an increase of almost ten times in 2031. This is a serious threat; a large mining area disturbs the karst landscape, thus affecting geomorphological conditions and eliminating the hydrological, biodiversity, climate change, and social functions of karst (He et al., 2021). Various forms of illegal mining in the Maros-Pangkep karst area require serious attention because several large mining areas for cement companies are officially recognized by the government and are located outside the Bulusaraung FMU. Apart from legal mines, there are illegal mines in the Maros-Pangkep karst area as well, for example, collecting limestone for building materials. Mining, which disregards the environment, negatively impacts karst sustainability.

Based on the land cover change analysis, it is necessary to anticipate a decrease in the forest cover area due to increased illegal mining and agricultural land. Mining in karst areas should be conducted within a specific limited area and set at a location with minimal negative environmental impacts. Legal mining has been accommodated; thus, the government should no longer open opportunities for other forms of illegal mining outside the permitted areas. Additionally, the deforestation and forest degradation of the karst Maros-Pangkep will impact the role of karst as a water reservoir, affecting the sustainability of springs from the karst area. In addition, post-mining land rehabilitation must be considered so that the succession process in disturbed areas may occur faster and benefit ecological and economic functions. A study on a karst area in China showed that successful reforestation in a karst area restored the flow of springs after they did not flow due to deforestation (Kang et al., 2020). Furthermore, agricultural practices already carried out in karst areas should apply soil conservation principles and prioritize mixed cropping patterns with tree species of non-timber forest products that have economic value. The ongoing forest rehabilitation implemented at the Bulusaraung FMU may foster land cover and increase biomass in the future. Apart from important research findings that we gained, future researches are needed using other model approaches with different parameters other than the CA-Markove Module and weighted each parameter in the study, hence it might reveal more interesting findings.

5. Conclusions

The secondary forests of the Maros-Pangkep karst ecosystem play a significant role as carbon sinks because they contain high biomass, carbon storage, and carbon sequestration. The study identified the karst plain contained 166.20 ton ha⁻¹ of biomass, 78.11 ton ha⁻¹ carbon stock, and 286.41 C ton ha⁻¹ carbon sequestration, respectively. Meanwhile, tower karst contains 160.65 ton ha⁻¹ biomass, 75.51 ton ha⁻¹ carbon stock, and 276.85 ton ha⁻¹ carbon sequestration. Unfortunately, the role of the secondary forest of the Maros-Pangkep karst hills as carbon storage is facing threats due to a tendency toward deforestation accompanied by an increase in mining and agricultural land. Forest encroachment and illegal limestone mining in karst area occurred intensively, which led to deforestation. The research findings suggest to control the opening of mining areas according to sustainable environmental principles, to encourage the implementation of post-mining land rehabilitation and to implement soil conservation rules by prioritizing agroforestry farming patterns and fostering forest rehabilitation.

6. Acknowledgement

The authors acknowledge the Asian Forest Cooperation Organization (AFoCO), the Center for Sustainable Forest Management Standardization Instruments, Ministry of Environment and Forestry of Republic of Indonesia, and Korea Forest Service (KFS) of Republic of Korea for development assistant and financial support of the research under a project collaboration “Innovative solutions for climate change and biodiversity landscape strategy to support SDGs in Indonesia” (project code: AFoCO/023/2021). We thank the Forest Management Unit of Bulusarung, South Sulawesi Provincial Forestry Service for their kind support in the implementation of project activities.

Appendix 1.

Two types of karst at the Maros

Figure appendix 1. Two types of karst at the Maros – Pangkep area. a. Plain karst; b. Tower karst.

Appendix 2.

T-test analysis results

T test carbon vegetation

Table

Independent Samples Test
		Levene’s Test for Equality of Variances		t-test for Equality of Means
		F	Sig.	t	df	Sig. (2-tailed)	Mean Difference	Std. Error Difference	95% Confidence Interval of the Difference
									Lower	Upper
Carbon	Equal variances assumed	.078	.782	.251	41	.803	2.94178	11.72973	−20.74688	26.63045
	Equal variances not assumed			.251	40.276	.803	2.94178	11.72201	−20.74423	26.62779

T test soil carbon

Table

Independent Samples Test
		Levene’s Test for Equality of Variances		t-test for Equality of Means
		F	Sig.	t	df	Sig. (2-tailed)	Mean Difference	Std. Error Difference	95% Confidence Interval of the Difference
									Lower	Upper
Soil_carbon	Equal variances assumed	.540	.466	−1.301	41	.200	−170.18546	130.80003	−434.34148	93.97056
	Equal variances not assumed			−1.279	35.845	.209	−170.18546	133.06800	−440.10051	99.72960

T test herb sequestration

Table

Independent Samples Test
		Levene’s Test for Equality of Variances		t-test for Equality of Means
		F	Sig.	t	df	Sig. (2-tailed)	Mean Difference	Std. Error Difference	95% Confidence Interval of the Difference
									Lower	Upper
Herb_sequestration	Equal variances assumed	2.561	.117	.880	41	.384	.36772	.41775	−.47595	1.21138
	Equal variances not assumed			.855	31.713	.399	.36772	.43019	−.50886	1.24430

T test sapling sequestration

Table

Independent Samples Test
		Levene’s Test for Equality of Variances		t-test for Equality of Means
		F	Sig.	t	df	Sig. (2-tailed)	Mean Difference	Std. Error Difference	95% Confidence Interval of the Difference
									Lower	Upper
Sapling_sequestration	Equal variances assumed	.432	.515	.248	41	.806	1.79915	7.26354	−12.86987	16.46818
	Equal variances not assumed			.252	40.689	.802	1.79915	7.14640	−12.63663	16.23494

T test pole sequestration

Table

Independent Samples Test
		Levene’s Test for Equality of Variances		t-test for Equality of Means
		F	Sig.	t	df	Sig. (2-tailed)	Mean Difference	Std. Error Difference	95% Confidence Interval of the Difference
									Lower	Upper
Pole_sequestration	Equal variances assumed	.406	.527	.109	41	.914	1.48216	13.62045	−26.02490	28.98922
	Equal variances not assumed			.107	35.892	.915	1.48216	13.85456	−26.61913	29.58344

T test tree sequestration

Table

Independent Samples Test
		Levene’s Test for Equality of Variances		t-test for Equality of Means
		F	Sig.	t	df	Sig. (2-tailed)	Mean Difference	Std. Error Difference	95% Confidence Interval of the Difference
									Lower	Upper
Tree_sequestration	Equal variances assumed	.602	.442	.133	41	.895	5.90914	44.56361	−84.08890	95.90717
	Equal variances not assumed			.134	40.983	.894	5.90914	44.17605	−83.30732	95.12560

T test Total sequestration

Table

Independent Samples Test
		Levene’s Test for Equality of Variances		t-test for Equality of Means
		F	Sig.	t	df	Sig. (2-tailed)	Mean Difference	Std. Error Difference	95% Confidence Interval of the Difference
									Lower	Upper
Total_sequestration	Equal variances assumed	.019	.895	.047	8	.964	3.82200	81.12738	−183.25808	190.90208
	Equal variances not assumed			.047	7.993	.964	3.82200	81.12738	−183.28612	190.93012