Optimization of Land Reuse Structure in Coal Mining Subsided Areas Considering Regional Economic Development: A Case Study in Pei County, China

Abstract

Land subsidence, which has caused large-scale settlement loss and farmland degradation, was regarded as the main constraint for land reclamation in the High Groundwater Coal Basins (HGCBs) in the eastern China plain. Both coal mining and agricultural production are important for regional development in this area. In general, the land reclamation direction in this area is greatly affected by the adequacy of filling materials and the land use demand of regional economic development. Taking seven coal mining subsided areas in Pei county, located in the eastern China plain, for example, this study proposed an integrated model (including the Limit Condition model, Logistic Regression model, Grey Linear Programming model and the conversion of land use and its effects at small regional extent (CLUE-S) model) to simulate and optimize the post-mining land use structure to meet the economic development needs of Pei county. Then, the post-mining land use structure under different scenarios, which were set based on the subsidence depth, were compared to explore the optimal collapse depth for separating the damaged land into the filling area and non-filling area. The landscape structure, ecological benefits, engineering quantity and reclaimed farmland area were used to compare the reclaimed land use structure of different scenarios. The results showed that the integrated model could efficiently simulate the reclaimed land use structure to meet the land demand for regional development. The optimal collapse depth for separating the damaged area into the filling area and non-filling area was 2.6 m. Currently, the reclaimed land use structure not only needs a low quantity of filling material, but also obtains a good landscape structure and elevated ecosystem services value. Furthermore, the reclaimed urban land was mainly distributed around Pei town, and the reclaimed farmland was mostly distributed in the area between Pei town and Weishan lake, which were consistent with the pattern of urbanization. The study provides valuable information for future land use plans for Pei county and will contribute to the methods of post-mining land reclamation in other HGCBs.

1. Introduction

It has been acknowledged that the mining industry supported national economic development, but also caused a vast amount of negative impacts [1], such as soil production capacity loss [2], land occupation [3], air and water pollution [4], geological disasters [5], biodiversity reduction [6], etc. China certainly could not avoid these problems. These problems were even worse in the High Groundwater Coal Basins (HGCBs), where the mining industry and agricultural production were equally important for local development (called “Overlap Regions”) [7]. Underground mining in HGCB areas has disturbed the surface severely and resulted in large-scale land subsidence, leading to much ground settlement loss and farmland degradation [8,9,10]. In addition, as overlap regions are coupled with large populations, land subsidence also aggravated the conflicts between land use and local social and economic development [11]. Land reclamation was commonly considered as an effective tool to address the conflict between mining activities and farmland protection [12]. Suitability evaluation or determining optimal land use is the key point for land reclamation [13]. However, the complex and multiple problems in HGCB regions required land reclamation projects to not only resolve the land degradation and loss problems but also to arrange residents’ living space and create a good land use ecosystem. Moreover, as land reclamation plays an important role in adjusting the land use structure and land use has entered the era of stocks in China [14], land rehabilitation is playing an increasingly important role in supporting national food security, agricultural modernization and the coordinated development of urban and rural regions [12,15,16]. Therefore, it is necessary to develop a system of post-mining land use to support the needs of production, living and ecology demands in HGCBs.

A plethora of methods have been developed for land use suitability evaluation. For example, Cheng et al. [13] evaluated the reclamation suitability of damaged mine land using the integrated index and the difference-product method. Palogos et al. [17] designed the optimal land use for reclamation based on evolutionary algorithms. Bangian et al. [18] demonstrated that the fuzzy decision-making method could efficiently optimize post-mining land use for pit areas. In addition, the multi-criterion decision-making methods and Analytic hierarchy process decision support system were used for determining the optimum manner of post-mining land use by Soltanmohammadi et al. [19] and Bascetin [20], respectively. However, most of these assessments or comparisons rarely avoid the subjectivity of evaluators. In addition, for an open mine, few studies designed the post-mining land use by considering the surrounding land use change during long-term mining activities. As for post-mining land use directions, agricultural land, forest land, and grass land are still the main land reclamation directions in China according to the Measures for the Implementation of the Regulation on Land Reclamation. However, agricultural land alone can never resolve the social, economic and environment problems in HGCB regions. In fact, various kinds of land use type, such as forestry, agriculture, commercial, recreation and residential, were recommended after mining activity ceased in the USA and Poland. Kazmierczak et al. [21] indicated that the most important aspect concerning the restoration of post-mining areas’ utility was determining the functions they should have after the mining activity ceased. Lima et al. [22] also stated that the land reclamation objective after mine closure should be site-specific. Therefore, it is reasonable to explore multifarious post-mining land use in HGCB regions.

In addition, unlike the damaged land in other areas, land subsidence is the dominating factor impacting the environmental and social issues in HGCBs [9], which makes the post-mining land use highly dependent on the adequacy of the filling materials. Without filling materials, subsidence areas with high collapse depths can only be converted to pond, whereas post-mining land use can be diverse with enough filling soil. As a result, in this study, we try to optimize the reclamation direction for coal mining subsidence land in HGCBs considering the adequacy of filling materials and the land use demand driven by socio-economic development. Seven coal mining subsidence areas in Northern Pei county were studied in the current research. According to the technical system created previously [23], this study has constructed an integrated model consisting of the Limiting method, Logistic Regression model, Grey Linear Programming model and CLUE-S (the Conversion of Land Use and Its Effects at Small Regional Extent) model, to explore post-mining land use structure. This structure was expected to be coordinated with the land use structure of Pei county and support its social and economic development. In addition, several scenarios of filling the damaged areas with different collapse depths were set up to explore the optimal reclamation filling depth. The results would provide an important basis for the optimization method for determining post-mining land use structure in HGCBs.

2. Materials and Methods

2.1. Study Area

Seven coalfields in the Northern Pei county, Jiangsu province, China (longitude 116.70° E–117.16° E, latitude 34.46° N–34.94° N), were selected as the study area, including Longdong coal mine, Yaoqiao coal mine, Xuzhuang coal mine, Kongzhuang coal mine, Zhangshuanglou coal mine, Longgu coal mine and Peicheng coal mine. These mine sites were located in the Huanghuai alluvial plain with high groundwater levels and dense villages and are typical overlap regions of coal mining and agricultural production. The soil layer in this area is thick and belongs to tidal soil. The climate is the north temperate with an average annual precipitation of 789.20 mm that is mainly concentrated in June and September. Therefore, it is prone to drought in spring and waterlogging in summer [24]. Since mining in 1977, coal mining has caused a large area of land subsidence, water accumulation and settlement damage, resulting in serious conflicts between human activity and land preservation. According to the mining plan of coal resources of Pei county, it was estimated that 8273.43 ha of land in the study area would subside in 2030 (Figure 1), which would be distributed across Longgu town, Yangtun town, Datun town, Lulou town, Anguo town and Pei town (Figure 1a). The predicted maximum collapse depth was 4.20 m (Figure 1b). As the depth of groundwater level is about 1.5 m, the overlapping land (e.g., farmland) would be readily changed into pond or wetland. Therefore, the land subsidence and subsequent damage to the local living environment in this area would seem to be irreversible.

2.2. Data Source

In this study, the land use maps of Pei county in 2006 and 2015 were used to train and validate the CLUE-S model. Referring to the Land Use Classification in China (GB/T 21010-2007) and research needs, the land use map was reclassified into nine categories, including cultivated land (CL), garden land (GL), forest land (FL), other agricultural land (OAL), urban land (UL), rural residential land (RRL), traffic and water conservancy land (TWL), other construction land (OCL) and water area and natural reserve land (WNL). Brief explanations of the different land use types were given by Chen and Zhou [25]. Eleven driving factors for land use change were selected, including four socio-economic factors, two natural driving factors and five location factors. The data of social and economic factors were obtained from the statistical yearbook of Pei county. The natural driving factors were collected from the digital elevation model of Pei county, with a spatial resolution of 30 m. The location factors were extracted from the land use map. The 13th Five-Year Plan outline for national economic and social development and the general land use plan (2006–2020) of Pei county were used to predict the land use quantity structure of Pei county in 2030 for meeting the demands of economic development. The mining engineering map, mining plan and collapse prediction parameters of each mining area were used to predict the land collapse distribution map of the study area in 2030 using a probabilistic integration method.

In this paper, ArcGIS 10.5 was used to reclassify the land use maps, extract the driving factors, analyze the balance of excavation and filling earth in the mining subsidence area and make the spatial distribution map. FRAGSTATS 4.1 was used to analyze the landscape pattern. SPSS 20.0 was used to implement the logistic model.

2.3. Method

2.3.1. Overview of the Integrated Model

Generally, the adequacy of filling materials determines the reutilization direction of reclaimed land in HGCBs. Therefore, this integrated model included four parts. Firstly, the mining subsidence area was divided into a filling area and non-filling area based on the sufficiency of the filling materials and the depth of subsidence. Secondly, the land use suitability evaluation was implemented in the non-filling and non-subsidence areas, and filling area, respectively. On one hand, the reclamation direction of the non-filling area was determined as other agricultural land, such as fishpond, using the Limit Condition method due to the deep accumulated water in the non-filling area. On the other hand, the Logistic Regression model was used to predict the suitability probability of the filling area and non-subsidence area for each land use type. Thirdly, the Grey Linear Programming model, which is a method of grey systems analysis for decision making under uncertainty [26], was used to predict the land use quantity structure of Pei county in 2030 to satisfy the land use demand of regional economic development. Lastly, the land use spatial structure of Pei county including the collapsed area under different scenarios was simulated by the CLUE-S model based on the land quantity structure demand and the land use suitability evaluation results. The optimal land reutilization quantity and spatial structure of the collapsed area were determined by comparing the ecological benefits, landscape structure, cultivated land area and reclamation engineering quantity of reclaimed land in different scenarios, which were the main concern indices in the current evaluation system of the land reclamation project in China. The simulation scenarios were designed by the depth of collapse. The technical framework is shown in Figure 2.

2.3.2. Land Use Suitability Evaluation

• Non-filling Area

The suitable reclamation direction of the non-filling area was determined as other agricultural land such as fishpond or wetland using the Limit Condition method, as the soil function was restricted by water accumulation in this area. Its area was set as x₀ (Equation (1))

$$x_0=A|h\ge h_m$$

(1)

where h is the collapse depth, h_m is the collapse depth that separates the damaged land into a non-filling area and filling area and A represents the area of the non-filling area.

• Filling Area and Non-Subsidence Area

The Logistic Regression model was used to evaluate the land use suitability (P_j,u) for each land use type considering the natural and socio-economic characteristics of each evaluation unit. The calculation formula (Equation (2)) was as follows:

$$Log\left(\frac{p_{j,u}}{1-p_{j,u}}\right)={\beta }_0+{\beta }_1{x}_{1,j}+{\beta }_2{x}_{2,j}+\cdots +{\beta }_n{x}_{n,j}$$

(2)

where P_j,u represents the probability of the jth evaluation unit for land use type u, x_n,j represents the value of the jth unit for the nth factor and β_n represents the regression coefficient of the nth factor.

The reliability of Logistic Regression model was assessed by the Relative Operating Characteristic (ROC) method, in which the area below the ROC curve was used to check the prediction accuracy. It is demonstrated that 0.7 < ROC < 0.9 represents an excellent prediction accuracy, while 0.5 < ROC < 0.7 indicates a low prediction accuracy. If the ROC value is bigger than 0.9, the prediction results are very reliable [27].

2.3.3. Optimization of Land Use Quantity Structure

There are many models used for land use quantity structure optimization, such as the Linear programming model [28], Markov chain [29], System dynamics model [30] and Neural network [31]. However, except for the Linear programming model, most methods predict land use structure through drawing upon past land use change practices. The socio-economic development of resource-based cities is highly affected by the mineral resources. With the consumption of mineral resources, the economic structures of resource-based cities will also be changed. Therefore, the Grey Linear Programming model, which is not based on the land use change history and permits certain gray component in optimization procedure, was proposed to predict the land demand of Pei county in 2030.

(1) Optimization Goal

Taking the maximization of land use economic benefits as the objective and the area of each land use type as the independent variable, the objective function (Equation (3)) was defined as follows:

$$f(x)=\max \left({\displaystyle \sum _{i=1}^3{c}_i\times x_i+c_4\times \left(x_0+x_4\right)+{\displaystyle \sum _{i=5}^9{c}_i\times x_i}}\right)$$

(3)

where f(x) expresses the land use economic benefits, 10,000 CNY/ha; x_i represents the area of the ith land use type and x₁–x₉ represent the CL, GL, FL, OAL, UL, RRL, TWL, OCL and WNL, respectively; and c_i represents the corresponding economic benefit coefficient of the ith land use type, 10,000 CNY/ha.

The economic benefit coefficients of different land use types were obtained from previous studies and the statistical yearbook of Pei county (

Table 1. The economic benefit coefficients of different land use types (10,000 CNY/ha).

Land-Use Types	CL	GL	FL	OAL	UL	RRL	TWL	OCL	WNL
Economic Benefit Coefficient	7.11	21.45	4.35	18.46	832.44	1	165.46	22.93	1

).

(2) Constraint Condition

Referring to previous studies, constraints were set up from seven aspects: the total land area constraint, land use planning constraint, land development intensity and conservation and intensive use constraint, labor resource constraint, land use ecological benefits and social benefits constraint, and model constraint. All of the constraints are shown in

Table 2. Expression and illustration of the limiting factors for land use quantity structure optimization in Pei county.

Constraints Type	Formula Expressions	Constraints Explanation
Total land area constraints	$x_1+x_2+\cdots +x_9+x_0=\mathrm{180,577.70}$	The total area of all kinds of land and filling area was fixed.
Land use planning constraints	$x_1+x_2+x_3+x_4+x_0\ge \mathrm{103,092.79}$ $x_4+x_0\ge \mathrm{14,771.49}$ $x_2\ge 4532.04$ $x_3\ge 1189.41$ $x_1\ge \mathrm{72,093.33}$ $762.87\le x_8\le 875.92$	The agricultural land should not less than the current area, and the garden land and forest land nearly keep steady according to the Land Use Plan of Pei county from 2006 to 2020; The cultivated land should not less than the current area in accordance with farmland protection policy; The upper and lower limits of other construction land were determined based on the 13th Five-Year Plan for coal industry development and the coal resource mining plan of Pei county.
Land development intensity and conservation and intensive use constraints	$x_5+x_6+x_7+x_8\le \mathrm{27,567.92}$ $x_5+x_6\le \mathrm{24,515.91}$ $m_1\le \frac{x_7}{x_5+x_8}\le m_2$ $x_5\le P_2\times a_2$	The total area of construction land should less than the current area complying with the principle of saving and intensive land use; The area of urban land and rural residential land should not more than the current area according to the regulation of the general land use plan of Pei county; The upper limit of urban land area was determined according to the standard for village and town planning (GB20188-1993); The area ratio of traffic and water conservancy land to the area of urban land and other construction land should not be less than the current ratio and should be higher than the ratio of adjacent developed cities.
Labor resource constraints	$b_1{x}_1+b_2{x}_2+b_3{x}_3+b_4\left(x_4+x_0\right)\le P_1$	The labor force for agricultural land should not higher than the rural population.
Land use ecological benefits constraints	$x_9\ge \mathrm{50,058.5}$	According to the ecological red line of Pei county, the area of water area and natural reserve land should not be less than 50,058.5 ha.
Land use social benefits constraints	$p_1\times a_1\le x_6\le \mathrm{18,598.99}$	The low limit area of rural residential land was determined by the standard for village and town planning (GB20188-1993), satisfying the living demands of people, and the upper limit was set as the current area, adapting the trend of urbanization.
Model constraints	$x_i\ge 0(i=0,1,2\cdots ,9)$	All variables in the model should be higher than zero.

Note: p₁, p₂ represent the rural and urban population in 2030, respectively, which were deduced by the natural population growth index method, and the natural population growth rate was set as 8‰ referring to the statistic book of Pei county. m₁ and m₂ represent the lower and upper limits of the ratio of traffic and water conservancy land area and the area of urban land and other construction land, which were set as 0.46 and 0.55, referring to this proportion of Xuzhou city, respectively. b₁, b₂, b₃ and b₄ represent the labor force coefficients of cultivated land, garden land, forest land, and other agricultural land, which were set as 2.66, 7.5, 7.5 and 3.61, respectively [32].

.

2.3.4. Optimization of Land Use Spatial Distribution

The spatial allocation module of the CLUE-S model was used to optimize the land use spatial distribution of Pei county and the collapsed area. The principle of this module includes three steps and can also refer the study of Peter H. Verburg [33].

(1)

Recognize the unit whose land use type can be changed, and the unchanged unit will be excluded before simulation for spatial allocation.

(2)

Calculate the total probability of the jth unit for the land use type u (Equation (4)).

$$TPRO{P}_{j,u}=p_{j,u}+ELA{S}_u+ITE{R}_u$$

(4)

where P_j,u refers the suitability probability of the jth evaluation unit for land use type u; and ELAS_u and ITER_u represent the elastic coefficient and the iterative variable of land use type u, respectively.

(3)

Spatial allocation. The total suitability probability of each unit for each land use type was calculated, and then the land use type with the bigger total probability was selected for a certain unit.

The simulation accuracy of the CLUE-S model was evaluated according to the Kappa coefficient value, which, if close to 1, indicated good similarity between the actual land use map and simulated land use map [34].

2.3.5. Scenario Simulation and Comparison

The sufficiency of the filling materials affects land reclamation in the mining subsidence area. Generally, the “Digging Deep to Fill Shallow” method obtained widespread availability [1]. If there are not enough external filling materials, the filling engineering will be carried out gradually from shallow to deep. The collapse depth, which separates the damaged land into the filling area and non-filling area, was 1.8 m when the excavation and filling earthwork were balanced in the case of no external filling materials through the cut and fill analysis. In the current study area, the depth of ground water is about 1.5 m and the depth of a fishpond is normally 2.5 m. Therefore, the optimal reclamation direction of a collapse unit is fishpond if the collapse depth reaches more than 4 m. As a result, the adequacy of external materials determines whether the area with a collapse depth between 1.8 m and 4.0 m can be filled. Considering the sufficiency of the filling materials, we designed eight scenarios based on different collapse depths (1.8 m, 2.0 m, 2.2 m, 2.4 m, 2.6 m, 2.8 m, 3.0 m and 4 m) separating the subsidence area into the filling area and non-filling area.

The ecological benefits, engineering quantities, landscape structures and cultivated land areas of the reclaimed land use under different scenarios were compared. The ecological benefits of the reclaimed land could be represented by the reclaimed land use ecosystem service value. The standard ecosystem service value was 2171.17 yuan (RMB)/ha [35], calculated according to Xie et al. [36] and Zhang et al. [37]. The engineering quantities were represented by the quantity of external earthwork. The landscape structure was represented by a series of landscape type indices and landscape diversity indices [38]. However, a number of optimal landscape type indices (NOLTI) for each scenario as an integrity index was compared, as no single landscape type index could represent the landscape structure totally. Briefly, the first step was to calculate and compare the landscape type index value of different reclaimed land use type in different scenarios. Then, when the landscape type index was the largest, the number of optimal landscape type indices for the corresponding scenario was counted as 1 (otherwise counted as 0). Finally, the number of optimal landscape type indices for each scenario was acquired. The area of reclaimed land was represented by the quantity of cultivated land after reclamation. All indices are shown in

Table 3. The index system for assessing the reclaimed land use structure of the subsided area.

Items	Index	Index Meaning
Ecological benefits	Ecosystem service value	The supply, regulation, support, culture and other service functions that the ecosystem can provide.
Landscape structure	Largest Patch Index (LPI)	The proportion of the largest patch area of a landscape type to the total landscape area.
	Proportion of Like Adjacencies (PLADJ)	The adjacency ratio between patches of a certain landscape type.
	Interspersion Juxtaposition Index (IJI)	The proximity between patches of one landscape type and patches of other landscape types.
	CONNECT	The adjacency between patches in a landscape type.
	DIVISION	The degree of patch separation in a landscape type.
	Aggregation (AI)	The spatial aggregation of patch in a landscape type.
	Shannon’s Diversity Index (SHDI)	The complexity of landscape structure composition and landscape heterogeneity.
	Shannon’s Evenness Index (SHEI)	Whether the landscape is dominated by one or a few dominant patch types.
Engineering works	External earthworks	The volume of external filling materials.
Cultivated land quantity	The area of reclaimed cultivated land	The area of reclaimed cultivated land in the subsidence area.

.

3. Results

3.1. Distribution of the Filling and Non-Filling Areas in the Mining Subsidence Area under Different Scenarios

The distribution of the filling and non-filling areas is shown in Figure 3. From scenario 1 to scenario 8, the area with a collapse depth of less than 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0 or 4.0 m was filled with exogenous materials, while the other area was left as a pool zone, respectively. It can be seen that the area of the non-filling area gradually decreased with the increase in excavation depth and external filling materials. This means that more and more subsided land can be reclaimed as other land categories rather than only as pool.

3.2. Land Use Quantity Structure of Pei County in 2030

The land use quantity structure of Pei county in 2030 predicted by the Grey Linear Programming model under different scenarios is shown in

Table 4. The land use quantity structure of Pei county in 2015 and 2030.

Year	CL	GL	FL	OAL	UL	RRL	TWL	OCL	WNL
2015	82,599.85	4532.04	1189.41	14,771.49	5024.55	18,598.99	2731.16	875.92	50,254.29
2030	72,093.33	4532.04	1189.41	26,386.94	8211.73	12,231.62	4998.21	875.92	50,058.50

. Due to the small area of the non-filling area, the land use quantity structure of Pei county in 2030 was the same under different scenarios. Compared with the land use quantity structure in 2015, the area of cultivated land decreased due to the small economic benefit coefficient of cultivated land. However, the area of cultivated land still met the requirement of cultivated land retention according to the general land use plan of Pei county. The increase in urban land and the decrease in rural residential land were consistent with the development trend of urbanization and the requirements of rural land renovation. The increase in traffic and water conservancy land complied with the requirements of economic development. The increase in other agricultural land also met the requirements of the general land use plan of Pei county.

3.3. Land Use Spatial Structure of the Subsidence Area

3.3.1. Validation of the CLUE-S Model

The Logistic Regression coefficients are shown in

Table 5. The Logistic Regression coefficients of different driving factors for each land use type in Pei county.

Driving Factors	Land Use Type
	CL	GL	FL	OAL	UL	RRL	TWL	OCL	WNL
Elevation	0.9410	1.0560	1.0490	0.9930	-	1.1260	-	1.0260	0.9700
Slope	1.0080	-	1.0300	1.0070	-	0.9900	1.0180	1.0240	0.9540
Population	1.0000	1.0000	1.0004	1.0000	-	-	1.0000	1.0000	1.0000
Population density	1.0018	-	0.9706	1.0018	-	-	1.0026	1.0031	0.9978
GDP per town	1.0000	1.0000	0.9999	1.0000	-	1.0000	1.0000	1.0000	1.0000
GDP per hectare	0.9994	1.0000	1.0096	0.9994	-	1.0000	0.9992	0.9990	1.0007
Distance to river	1.0003	1.0002	1.0006	1.0001	-	1.0003	0.9999	0.9997	0.9982
Distance to first level road	0.9996	0.9999	1.0002	1.0002	-	0.9995	0.9943	-	1.0008
Distance to rural road	0.9983	0.9991	0.9988	0.9989	-	0.9997	1.0009	1.0009	1.0030
Distance to town	1.0000	0.9998	0.9996	1.0000	0.7210	-		-	1.0002
Distance to mining and industry land	1.0000	1.0001	1.0007	0.9998	-	0.9999	1.0000	0.9964	1.0000
Constant	15.29	0.01	0.00	0.23	32,462.14	0.01	0.07	0.03	0.08
ROC	0.78	0.74	0.95	0.69	0.99	0.76	0.91	0.95	0.96

. Except for OAL, the Receiver Operating Characteristic (ROC) curve values corresponding to each land use type were above 0.7, indicating that the selected driving factors could be used to simulate the land use change of Pei county [39]. By using a trial and error method, the elastic (ELAS) coefficients of CL, GL, FL, OAL, UL, RRL, TWL, OCL and WNL in the CLUE-S model were determined as 0.6, 0.4, 0.5, 0.35, 0.1, 0.45, 0.3, 0.2 and 0.6, respectively. Besides, the kappa coefficient between the simulated land use map in 2015 and the actual land use map in 2015 was 0.93, indicating that the corrected CLUE-S model can be used for scenario simulation analysis.

3.3.2. Comparison of Different Scenarios

The ecological benefits, engineering quantities, landscape structure and cultivated land areas of the reclaimed land use under different scenarios were shown in

Table 6. A comparison of the reclaimed land use structure of the coal mining subsided area in different scenarios.

Scenario	SHDI	SHEI	NOLTI	Reclaimed Cultivated Land/ha	External Earthwork/106 m3	Ecosystem Service Value/108 CNY
Scenario 1	1.3794	0.6278	1	995.40	0.00	11.45
Scenario 2	1.3877	0.6315	4	1096.20	11.46	11.55
Scenario 3	1.4013	0.6378	5	1241.64	24.71	11.48
Scenario 4	1.4048	0.6394	7	1296.72	34.56	11.46
Scenario 5	1.4140	0.6435	12	1401.84	42.21	11.42
Scenario 6	1.4127	0.6429	6	1369.08	47.41	11.44
Scenario 7	1.4131	0.6431	10	1337.04	50.05	11.46
Scenario 8	1.4113	0.6423	7	1320.48	52.64	11.43
Current in 2015	1.3667	0.6220	-	-	-	10.99

Note: SHDI: Shannon’s Diversity Index; SHEI: Shannon’s Evenness Index; NOLTI: Number of optimal landscape type metrics.

. Scenario 5 had the largest NOLTI, indicating that patches connectivity and cohesion in this scenario were better than those in other scenarios. Besides, scenario 5 also obtained the highest Shannon’s Diversity Index, Shannon’s Evenness Index and reclaimed cultivated area among all scenarios. The ecosystem service value was the lowest in scenario 5, but it was still higher than the value of the study area before reclamation. In fact, the ecosystem values of eight scenarios were all higher than their original values. This result was consistent with a previous study that reported the impacts of land use change on ecosystem services in HGCBs [7]. It mainly resulted from the increase in water area due to the subsidence.

Therefore, scenario 5 obtained the optimal reclaimed land use quantity and spatial structure. The reclaimed land use spatial distribution map is shown in Figure 4. The areas of reclaimed CL, GL, FL, OAL, UL, RRL, TWL, OCL and WNL were 1401.84, 40.32, 2.88, 1793.32, 184.32, 296.64, 27.72, 433.08 and 3763.96 ha, respectively. Among them, the reclaimed UL was mainly distributed in Pei town, which was in line with the trend of land expansion in the process of urbanization. In addition, other agricultural land was mainly near to Weishan Lake, and CL was mainly distributed between Weishan Lake and Pei town. Lastly, some other construction land was reserved, following the coal mining plan of Pei county.

4. Discussion

4.1. Simulation Accuracy of the Combined Model

In this study, ROC values and the kappa coefficient were used to test the prediction accuracy of the coupled model. The ROC values of CL, GL, FL, OAL, UL, RRL, TWL, OCL and WNL were 0.78, 0.74, 0.95, 0.69, 0.99, 0.76, 0.91, 0.95 and 0.96, respectively. It meant that the Logistical regression model achieved high prediction accuracy for FL, UL, TWL, OCL and WNL and satisfactory results for CL, GL and RRL. However, the ROC value for OAL was a little bit lower than 0.7, indicating a low prediction accuracy. This might be mainly attributed to the complex constitution of OAL, which was not recommended for Logistic Regression according to the guideline of the CLUE-S model [39]. Besides, the autocorrelation commonly has an impact on the ROC value. For example, Jiang et al. [40] reported that the ROC value for bare land was 0.63 based on Logistic Regression model, but increased to 0.833 based on the Auto-Logistic Regression model. Therefore, the ROC value of OAL may be elevated if Auto-Logistic Regression is used.

The kappa coefficient of the simulated land use map and actual land use map of Pei county in 2015 was 0.93, suggesting that the calibrated model could effectively capture the evolution trend of land use change in Pei county. This value was a little bigger compared with the kappa coefficient reported by other studies using the CLUE-S model [41,42,43,44]. However, Zhang et al. [45,46] found that the kappa coefficient value was around 0.90 in their two studies about land use change in Pinggu and Mentougou districts. Through comparing the characteristics of the above studied areas with those of our study area, the small scale and plain terrain of our study area may facilitate the relatively high prediction accuracy of the CLUE-S model. In the following research, the effect of the scale and terrain characteristic of a study area on the simulation accuracy of the CLUE-S model could be investigated.

4.2. The Collapse Depth Separating the Damaged Land into the Filling and Non-Filling Areas

Land subsidence usually causes water accumulation, sloped surfaces and fissures, which then affect soil density, soil water content and soil hydraulic conductivity [47]. In addition, land collapse often forms saline-alkali soil and saltwater lakes, especially when the land is close to the water area [47]. When the collapse depth is higher than the maximum groundwater depth, damaged land will become seasonally or perennially waterlogged [11]. In this study area, the maximum groundwater depth is about 1.5 m, indicating that crops cannot survive when the depth of the collapse is greater than 1.5 m. This study determined the optimum filling area through comparing the ecological benefits, engineering quantities and landscape structures of different scenarios. Subsidence areas with collapse depths higher than 2.6 m were recommended as fishpond, water reservoirs or other agricultural land. This result was consistent with an exploring policy which aimed to make a criterion to assess if the land use type on the cadaster could be changed. This policy serves to permit the alteration of land types in the third national land survey for those damaged by long-term environmental changes and mining activities. In fact, according to field survey, subsided areas with collapse depths higher than 3 m were kept as water areas directly in engineering design. Therefore, this study provides a basis for local governments to apply for land alteration in areas with predicted subsidence depths higher than 2.6 m and a point of reference for making land use type changes in HGCBs.

4.3. Implications for Further Research

Land use suitability assessment for the coal mining subsidence area and non-subsidence area was the key of this technical framework, and then the results were used for spatial allocation. In this study, subsidence was regarded as the only limiting factor for land reclamation, and the subsidence area was divided into the filling area and non-filling area according to the depth of collapse. Then, the land use suitability for the non-filling area and non-subsidence and filling area were evaluated by the Limit Condition method and Logistic Regression method, respectively. Finally, the combined land use suitability results were used for spatial allocation by the CLUE-S model. However, the soil fertility difference between the filling area and non-subsidence area was not included in the land use suitability evaluation system. The spatial suitability competition was the core of the CLUE-S model to realize the spatial allocation function [33]. As a result, further research on the integral land use suitability assessment are needed to develop a tightly coupled model for achieving the optimization of land reclamation structure.

5. Conclusions

In this study, we proposed an integrated model to simulate the reclaimed land use structure and determined the optimal collapse depth for separating the filling area and non-filling area. By comparing the reclaimed land use structure under different scenarios, the optimal land reuse structure and the optimal collapse depth were obtained, which not only satisfy the economic development of Pei county but also provide good landscape structure. Specific conclusions are listed as follows:

(1)

The proposed integrated model could optimize the post-mining land use structure to satisfy the land demand of regional economic development and address the conflicts between mining activity and land use in Pei county. The application and suitability of the method should be further verified in other post-mining land use optimizations in HGCBs.

(2)

In Pei county, keeping areas where the collapse depth was higher than 2.6 m as fishpond or other agricultural land and filling the remaining area not only reduced the cost of the filling materials but also achieved good landscape structure and improved ecosystem service value.

(3)

The reclaimed areas of cultivated land, garden land, forest land, other agricultural land, urban land, rural residential land, traffic and water conservancy land, other construction land and water area and natural reserve land in 2030 were 1401.84, 40.32, 2.88, 1795.32, 184.32, 296.64, 27.72, 433.08 and 3765.96 ha, respectively. Among those studied areas, the subsidence area near Weishan lake was mostly reused for other agricultural land, while the reutilized urban land was mainly distributed around Pei town, consistent with the trend of urbanization. The cultivated land was mostly distributed in the area between Pei town and Weishan lake.