Integrating System Dynamics, Land Change Models, and Machine Learning to Simulate and Predict Ecosystem Carbon Sequestration Under RCP-SSP Scenarios: Fusing Land and Climate Changes

Abstract

Understanding the impacts of land use and vegetation carbon sequestration under varying climate scenarios is essential for optimizing regional ecosystem services and shaping sustainable socioeconomic policies. This study presents a novel research framework that integrates a system dynamics (SD) model, a patch generation land use simulation (PLUS) model, and the random forest algorithm, coupled with SSP-RCP scenarios from Coupled Model Intercomparison Project Phase 6 (CMIP6), to simulate future vegetation net primary production (NPP). A case study in Hubei Province, central China, demonstrates the framework’s effectiveness in elucidating the interactions between land use change, climate change, topography, and vegetation conditions on carbon sequestration. The integration of SSP-RCP scenarios provides a clear understanding of how different climate conditions influence regional carbon sinks, offering valuable scientific insights for regional carbon neutrality and sustainable development policymaking. The simulation results for Hubei Province across the years 2030, 2040, 2050, and 2060, under three pathways—SSP1-1.9, SSP2-4.5, and SSP5-8.5—reveal that SSP1-1.9 leads to the highest carbon sequestration, while SSP5-8.5 results in the lowest. The annual total carbon sink ranges from 115.99 TgC to 117.59 TgC, with trends varying across scenarios, underscoring the significant impact of policy choices on local ecosystems. The findings suggest that under low-carbon emission scenarios, there is greater potential for NPP growth, making carbon neutrality goals more achievable.

1. Introduction

Global climate change and environmental pollution present significant threats to the stability and sustainable development of human society. The rise in atmospheric greenhouse gases is a primary driver of these challenges [1,2]. Reducing carbon emissions to achieve carbon neutrality is an effective strategy for mitigating climate change [3]. As global sustainable development strategies gain traction, more countries are actively engaging in carbon reduction efforts, each crafting its own path toward carbon neutrality [4,5]. From the perspective of carbon neutrality, future global carbon emissions must be less than or equal to carbon sinks [6]. Carbon sinks serve as a critical benchmark for developing carbon reduction strategies, underscoring the importance of scientifically predicting and simulating regional carbon sink capacities when formulating carbon neutrality policies [7]. However, due to the current lack of appropriate tools and methodologies, many local governments often estimate future carbon sequestration levels based on experience rather than quantitative scientific analysis. For instance, some local governments in China rely on estimations to determine forestry carbon sequestration when developing policies to address climate change [8]. There is a pressing need for a scientific approach to accurately predict and simulate future regional carbon sink patterns and totals, providing a foundation for low-carbon development policies at national, regional, and urban levels [9]. Consequently, predicting and simulating future carbon sinks has become an increasing focus for researchers across various disciplines globally.

Current research on plant carbon sequestration primarily uses NPP as a key indicator. NPP represents the total energy fixed by plants through photosynthesis per unit area, minus the energy expended in plant respiration [10]. Research on NPP generally focuses on two areas: first, assessing the current status and characteristics of NPP through methods such as field surveys [11,12], eddy covariance, and model estimation [13], which have yielded insights into carbon sequestration trends, spatial distribution, and driving factors [14,15,16]; and second, predicting future NPP trends and spatial distribution, an area that remains less explored.

Predicting and simulating future plant carbon sequestration typically rely on various models, including linear regression models, carbon emission coefficient models, and the Thornthwaite Memorial model. Linear regression models establish a relationship between NPP and time, predicting future changes through trend extrapolation [17]. However, these models do not account for the impacts of climate, social, and land use and cover change (LUCC), potentially leading to significant prediction errors. Additionally, they are limited to quantitative predictions without simulating spatial distribution. Carbon emission coefficient models, such as the widely used InVEST [18,19], estimate carbon sequestration by calculating the total carbon emissions (positive) or absorption (negative) based on land use types and their corresponding emission coefficients. By multiplying the area of different land uses by these coefficients, the model assesses the overall carbon balance [20]. While effective in reflecting the impact of LUCC on carbon budgets, this method does not consider climate effects or the spatial heterogeneity of carbon sequestration [21]. The Thornthwaite Memorial model, used for evaluating climate and water resources, estimates NPP based solely on climate data, utilizing parameters such as potential evapotranspiration (PET), effective precipitation, humidity index (HI), and vegetation adjustment factors (Fc) [22]. Although this model provides total NPP predictions, it cannot simulate spatial distribution and overlooks the impact of LUCC on carbon sequestration.

Research indicates that plant carbon sink dynamics are influenced by multiple factors, including land use types, climate, topography, and plant growth conditions [23,24]. For example, different land use types exhibit varying carbon sink capacities; forests and grasslands generally have higher carbon sequestration potential, whereas construction and unused lands may have lower capacities or even act as carbon sources [25,26]. Climate change directly affects plant growth, distribution, and species composition, thereby influencing the carbon absorption capacity of terrestrial ecosystems. Different plant species have specific environmental requirements, and climate shifts can alter their habitats, potentially leading to species extinction or migration [27,28]. For instance, some low-altitude plants may migrate to higher altitudes as temperatures rise. Topographic features also significantly impact plant growth and ecosystem functions; steep slopes with rapid water drainage and poor soil fertility may limit plant growth and photosynthesis efficiency, resulting in lower NPP. By contrast, flat, low-lying areas with better water storage and nutrient accumulation tend to support more vigorous vegetation and higher NPP. Thus, plant carbon sinks are not determined by a single factor but are the result of complex interactions among LUCC, climate variations, topography, and plant growth conditions. Accurately predicting future plant carbon sinks requires a comprehensive approach that integrates these elements, a gap that existing research has yet to fully address [29].

As computer science advances, the application of machine learning algorithms in ecological research has become increasingly prevalent. Machine learning, grounded in statistical principles and algorithms, enables computer systems to learn from large datasets, automatically recognizing patterns and extracting valuable insights for predictions or decisions [30,31]. Research shows that machine learning can uncover latent associations among factors in complex ecological systems and capture nonlinear relationships, providing deeper insights into how various environmental factors influence vegetation growth [32,33]. Applying machine learning to study the combined effects of land use, climate conditions, topography, and plant characteristics on vegetation NPP not only enhances our understanding of ecosystem dynamics but also improves the accuracy of NPP predictions. Therefore, leveraging machine learning to analyze changes in vegetation NPP offers new perspectives and solutions for addressing the challenges posed by global climate change and land use transformations.

In this study, we developed a research framework that integrates a system dynamics (SD) model, a patch generation land use simulation (PLUS) model, and a random forest model, coupled with SSP-RCP scenarios proposed by CMIP6, to simulate future vegetation NPP. This framework was applied in a rapidly developing province in central China to guide the formulation of sustainable development policies for the region. Additionally, three development scenarios—SSP1-1.9, SSP2-4.5, and SSP5-8.5—were designed to explore the impact of different climate and development pathways on regional carbon sinks.

2. Materials and Methods

2.1. Study Area

Hubei Province, located in central China, is traversed by the Yangtze River and spans an area of 18,590 square kilometers. As of 2023, it had an estimated population of approximately 58.44 million, with urban residents comprising 64.67% of the total. Geographically, Hubei transitions from China’s second to third terraces, with mountainous terrain dominating the west, north, and east, while the central and southern regions are characterized by the flat expanse of the Jianghan Plain (see Figure 1). The province is rich in vegetation, featuring both deciduous and evergreen broadleaf forests, and serves as a transitional zone between the botanical regions of eastern and western China. Hubei is also home to Wuhan, the largest city in central China, a major grain-producing area, and a crucial ecological barrier. The province’s diverse topography and ecosystems make it an ideal case study for generating insights applicable to similar regions worldwide. Additionally, research in this area can help validate the generalizability of ecological models across different natural environments.

2.2. Data Sources

(1) NPP Data: The NPP data were derived from the MOD17A3 data product, which aggregates annual vegetation NPP by summing all 8-day Net Photosynthesis (PSN) products (MOD17A2H).

(2) Land Use Data: Land use data were sourced from the Chinese Academy of Sciences; for this study, these were reclassified into six major classes: cropland, grassland, water bodies, forest land, built-up land, and unused land.

(3) Historical Climate Data: Historical climate data, such as rainfall and temperature, were sourced from the China Meteorological Data Service Center (https://data.cma.cn/, accessed on 11 November 2023). Future climate data came from CMIP, specifically including three paths—SSP1-1.9, SSP2-4.5, and SSP5-8.5—of the CAS-ESM2-0 model (https://data.tpdc.ac.cn/zh-hans/data/, accessed on 11 November 2023).

(4) Digital Elevation Model (DEM) Data: The digital elevation model data were sourced from NASA’s Shuttle Radar Topography Mission (SRTM) and provided by the Geographic Spatial Data Cloud (https://www.gscloud.cn/, accessed on11 November 2023).

(5) Point of Interest (POI) Data: POI data were obtained from Baidu Maps, representing the spatial distribution of various public infrastructures in 2010, such as shopping centers, restaurants, and schools.

(6) Other Statistical Data: Additional data, including demographics and economic indicators, were sourced from the Statistical Yearbook published by the Hubei Provincial Government and the “Hubei Province Land Spatial Overall Planning 2021–2035”, which outlines protected areas such as basic farmland and significant ecological conservation zones.

2.3. Research Methods

2.3.1. General Procedure

The variation in NPP is influenced by a combination of factors such as land use, climate change, slope, and elevation, and these factors do not have a linear relationship with NPP. Therefore, when simulating future NPP, it is essential to consider the effects of these factors and incorporate future land use and climate projections into the model. In this study, we first predicted the spatial and quantitative changes in land use, and based on these predictions, we simulated the future distribution of land use. Subsequently, we integrated future climate data to ultimately model the future NPP.

The NPP prediction model proposed in this study consists of several integrated modules, as depicted in Figure 2. The process begins with the construction of a System Dynamics (SD) model, which incorporates land use data, climate data, and socioeconomic development data, along with the Shared Socioeconomic Pathways (SSPs) provided by CMIP6. This model is used to forecast the total area of various land types under different future scenarios (Section 3.1). Next, a patch generation land use simulation (PLUS) model is trained using POI data, road and infrastructure data, and land use data. The LUCC probabilities are estimated using a random forest algorithm, with key model parameters optimized through a genetic algorithm (Section 3.2). The predictions from Section 3.1 are then used to simulate the spatial distribution of land use under different future scenarios. Following this, an autonomous Latin hypercube sampling model is developed to achieve spatially random and uniform sampling, integrating autonomous sampling with Latin hypercube sampling to prevent spatially concentrated sampling. Finally, a random forest (RF) model is trained using historical data on land use, elevation, climate, and NPP. After ensuring the model meets accuracy requirements, it is employed to simulate future regional NPP distribution, utilizing the LUCC simulation results from Section 3.2 and the SSP data (Section 3.3). Each iteration of model training includes accuracy checks to ensure the reliability and precision of the experimental results. Three different SSPs were used in this study, particularly the CAS-ESM2-0 model using the SSP1-1.9, SSP2-4.5, and SSP5-8.5 pathways; SSP1-1.9 represents a sustainable development path with a strong emphasis on ecological conservation and significant shifts in energy and industrial structures toward low-carbon development. SSP5-8.5, on the other hand, represents a high-emission development path that prioritizes the continued large-scale use of fossil fuels. SSP2-4.5 serves as a medium-carbon pathway, balancing elements of both extremes.

2.3.2. Construction and Simulation of Land Use SD Model

The core idea of SD is to analyze the interactions among various components within a system to uncover overall behavior patterns and underlying mechanisms [34,35]. SD is widely used to analyze complex systems such as ecosystems, economies, and social systems. The process of change in land use quantity and structure is a complex phenomenon influenced by multiple factors, including natural geography, sociocultural dynamics, economic development, and climate change [36,37]. In this study, four subsystems—economic, climate, population, and land—were constructed. Functional relationships were established based on the interactions among factors within each subsystem and between subsystems. By inputting data on the economy, climate, and population, the model predicts the structure of land use quantities (see Figure 3). The model includes 33 variables in total, comprising 8 state variables (SVs), 4 rate variables (RVs), 4 data variables (DVs), and 17 auxiliary variables (AVs) (see

Table 1. Evolution of land use quantity structure SD model main variables and variable types.

	Subsystem	Variable Name	Units	Variable Types
1	Population subsystem	Population	10,000 people	SV
2		Urban population	10,000 people	AV
3		Rural population	10,000 people	AV
4		Population change	10,000 people/year	AV
5		Urban population rate	%	RV
6		Population change rate	%	DV
7	Economic subsystem	GDP	One hundred million dollars	SV
8		Agricultural production	One hundred million dollars	AV
9		industrial production	One hundred million dollars	AV
10		Tertiary industry	One hundred million dollars	AV
11		Fixed investments	One hundred million dollars	AV
12		Residential investment	One hundred million dollars	AV
13		Investment in secondary and tertiary industries	One hundred million dollars	AV
14		GDP change rate	%	DV
15	Climatic subsystem	Annual temperature	°C	SV
16		Annual temperature change	°C/year	AV
17		Annual temperature change rate	%	DV
18		Annual precipitation	mm	SV
19		Annual precipitation change	mm/year	AV
20		Annual precipitation change rate	%	DV
21	Land subsystem	Unused land	km2	AV
22		Proportion Unused land rate	%	RV
23		Grassland	km2	AV
24		Woodland	km2	SV
25		Woodland change	km2/year	RV
26		Water	km2	SV
27		Water change	km2/year	RV
28		Construction land	km2	SV
29		Construction land change	km2/year	AV
30		Demand for urban construction land	km2	AV
31		Village construction land demand	km2	AV
32		Cultivated land	km2	SV
33		Cultivated land change	km2/year	AV

).

Economic subsystem: This subsystem focuses on GDP as its primary SV, which is significantly influenced by the GDP growth rate. The GDP growth rate can be derived from historical data for modeling and validation, while future projections are informed by relevant government documents and literature. From GDP, other key variables such as fixed asset investment, agricultural output, and industrial output are derived. Fixed asset investment is further subdivided into residential and secondary/tertiary industry investments, which influence the demand for urban and rural construction land, ultimately determining the total construction land demand. Industrial output, supported by urban construction land, affects its scale, while agricultural output is closely correlated with the area of arable land, which directly influences it.

Population subsystem: This subsystem centers on population as its primary SV, with population growth rate and urbanization rate as key RVs. The total population is divided into rural and urban populations, both of which impact the demand for urban and rural construction land. The urbanization rate specifically influences the urban population, which in turn affects the demand for construction land.

Climate subsystem: This subsystem incorporates rainfall and temperature as its primary SV, along with corresponding RVs. Historical data serve as the basis for these RVs, while future projections are derived from specific CMIP6 datasets. This subsystem directly influences environmental conditions for vegetation, thereby affecting changes in grasslands, water bodies, and forests.

Land subsystem: The land subsystem is the core of the SD model, encompassing six land use types: grasslands, forests, arable land, water bodies, construction land, and unused land. Major land use types such as arable land, construction land, water bodies, and forests serve as SVs, each supported by corresponding RVs. Minor land use types like unused land and grasslands function as AVs. The SVs are interconnected with AVs from other subsystems, reflecting their influence. Unused land has remained relatively stable over the years, minimally impacted by AVs from other subsystems, thus maintaining a consistent proportion over time. Grasslands act as a regulator of total land use within the system, calculated as the difference between the total area and the sum of the other five land use types.

The variable parameters were primarily determined through statistical methods and rate equations. Functional relationships between variables were established by fitting historical data to various functions, with the function yielding the highest R² value selected as the best representation of these relationships. Furthermore, model parameters were adjusted to account for the specific conditions in Hubei Province.

2.3.3. The PLUS Model

The PLUS (V2.0) framework consists of two primary modules: the first module utilizes the random forest algorithm within the Land Use Expansion Analysis Strategy (LEAS) framework, while the second module employs a multi-class random patch seed cellular automaton model (CARS) [38,39]. The core principle of this model is to simulate the influence of human activities and natural processes on land use, thereby predicting future land use patterns, with a particular emphasis on the patch generation algorithm. This algorithm divides the land into a series of small patches, each representing different land use types, such as cropland, forest, and urban areas, to simulate the LUCC process. Accuracy validation is primarily conducted by computing the kappa coefficient and the Figure of Merit (FoM) between simulated and observed data.

In the simulation process, we applied restrictions to the land-use type transition rules in the model based on different SSPs and local policies. For example, in the SSP5-8.5 scenario, no additional restrictions were imposed, allowing for mutual conversions between all land types. In the SSP2-4.5 scenario, locally designated ecological protection lines and permanent croplands were prohibited from converting to other land types. The SSP1-1.9 scenario, building upon SSP2-4.5, further restricted the conversion of ecologically important areas to other land types.

2.3.4. Sampling

The datasets used in this study are raster data. Given the complex and diverse terrain and the distinct ecological variations across regions in Hubei Province, the spatial location of sampling significantly impacts the training outcomes. To address this, the study integrated Latin hypercube sampling (LHS) with the bootstrap method [40,41,42].

The process begins by determining the number of samples to be extracted from the raster data, with each feature corresponding to a specific attribute of the raster. The value range of each feature is then divided into equal intervals. Based on the number of sample points, a Latin hypercube table is constructed—a two-dimensional matrix in which the number of rows and columns corresponds to the number of sample points. In each row and column, an interval is randomly selected, and a random value within that interval is generated as the sampling point. Each interval is chosen only once across rows and columns, ensuring that each row yields a unique sample. Finally, for each row, a sampling point is generated based on the corresponding interval and random value from the Latin hypercube table (see Figure 4).

Since Latin hypercube sampling selects only one point per region at a time, and this point may correspond to a “NO data” attribute, the sampling process incorporates the bootstrap method. Initially, the desired number of samples is set, and if a “NO data” value is encountered, it is excluded from the total count, and sampling is repeated until the required number of valid samples is obtained.

2.3.5. Random Forest Model for NPP Prediction

The random forest algorithm builds each decision tree by employing random sampling and feature selection, then integrates the predictions through a process of voting or averaging (see Figure 5). This method ensures high accuracy, robustness, and interpretability [43,44]. In this study, the decision trees are constructed using the CART algorithm, which determines the optimal splitting criterion for nodes based on the calculation of each attribute’s Gini coefficient [45]. The calculation process is as follows:

$$G\mathrm{ini}\left(D, \mathrm{a}\right)={\displaystyle \sum _{k=1}^{\left|y\right|}{\displaystyle \sum _{k^{\prime }\ne k}{p}_k}}{p}_{k^{\prime }}=1-{\displaystyle \sum _{k=1}^{\left|y\right|}{p}_k^2}$$

(1)

Here, $G\mathrm{ini}\left(D\right)$ represents the Gini index of node D, $k$ denotes the number of categories in feature $a$, and $P \left(i\right)$ denotes the proportion (probability) of samples belonging to the $i$ category within feature $a$. A lower Gini index indicates higher node purity, suggesting that the samples within the node are more likely to belong to the same category. Conversely, a higher Gini index reflects a more diverse sample distribution and lower purity.

The Gini coefficient is defined as follows

$$Gini\_index(D,a)={\displaystyle \sum _{v=1}^v\frac{\left|D^V\right|}{\left|D\right|}}Gini(D^v)$$

(2)

In this context, $Gini\_index(D, \mathrm{a})$ represents the Gini coefficient of feature a in sample $D$, $Gini(D)$ represents the Gini value of sample $D$, and $V$ represents the possible values within sample $D$. When the $Gini\_index(D, \mathrm{a})$ value is smaller, the purity efficiency of the node is higher, indicating that the attribute is more suitable to be used as the root node or closer to the root node.

Most of the feature data used for predicting vegetation NPP are continuous variables (see

Table 2. Main characteristic data of the random forest model.

Subsystem	Variable Name	Units
Climatic	Mean annual temperature	°C
	Mean annual precipitation	ml
	Maximum surface temperature	°C
	Minimum surface temperature	°C
	Mean annual wind speed	m/s
	Average annual relative humidity	%
Landform	Altitude	m
	Slope	°
Location	Latitude	°
	Longitude	°
LUCC	LUCC
Last NPP	Last NPP	gC/(m2 × a)

). However, decision trees in random forests are structured as binary splits. Therefore, before a feature can be used for classification, it must be converted into a binary variable. The process is as follows:

Suppose the feature sample set $D$ has a total of a CV. Within this feature sample set, there will be n (n ≤ a) distinct values. These values are first sorted in ascending order, denoted as {a1, a2, a3, …, an}. Based on split point $t$, $D$ can be divided into subsets, and the Gini values for each subset are calculated using formula 6. Weighted sums are then used to derive the Gini value of $D$. The smallest Gini value for $t$ is selected as the optimal split point for further calculations.

At the leaf nodes, the primary information consists of vegetation NPP. Since NPP is a CV, the average NPP of the leaf nodes after classification is stored as the final data used in constructing the random forest.

By training multiple decision trees that collectively form a forest, predicting the future NPP of vegetation becomes a regression problem. The final regression result is obtained by averaging the outcomes of all the decision trees. The formula is as follows:

$$\overline{h}\left(x\right)={\displaystyle \frac{1}{T}}{\displaystyle \sum _{T=1}^T\left\{h\left(x,{\theta }_t\right)\right\}}$$

(3)

In this context, $\overline{h}\left(x\right)$ represents the final classification result; T denotes the number of decision trees; $h\left(x,{\theta }_t\right)$ stands for the classification result of a single decision tree; and ${\theta }_t$ denotes mutually independent, randomly sampled variables.

Finally, the feature data that has little influence on the model is filtered out, and the model parameters are optimized.

3. Results

3.1. Prediction of Land Use Structure in Different Scenarios

3.1.1. Evaluation of the SD Model

A comparative analysis was conducted between the land use quantity structures predicted by the SD model for the years 2000, 2005, 2010, 2015, and 2020, and corresponding LUCC observational data were provided by the Chinese Academy of Sciences. The results demonstrate that for cultivated land, construction land, water bodies, and forested land, the trends predicted by the SD model closely align with the observational data. Although unused land showed some deviation due to its smaller base, the maximum relative error was 3.86% (

Table 3. Statistical table of errors between simulated and observed values of land use quantity in the SD Model (unit: km²).

Year		2000	2005	2010	2015	2020
Cultivated land	Deviation area	0.00	263.18	−467.5	101.74	−183.82
	Relative error	0.00%	0.39%	−0.70%	0.15%	−0.28%
Woodland	Deviation area	0.00	309.27	235.03	−123.16	−382.08
	Relative error	0.00%	0.33%	0.25%	−0.13%	−0.42%
Grassland	Deviation area	0.00	28.61	581.65	153.88	397.46
	Relative error	0.00%	0.41%	8.38%	2.23%	5.72%
Water	Deviation area	0.00	61.65	−229.89	122.03	−54.61
	Relative error	0.00%	0.47%	−1.86%	0.99%	−0.41%
Construction land	Deviation area	0.00	−251.11	−115.7	−245.65	229.19
	Relative error	0.00%	−4.67%	−1.70%	−3.08%	2.95%
Unused land	Deviation area	0.00	−13.06	−6.55	−9.82	−2.77
	Deviation area	0.00%	−3.86%	−1.70%	−2.57%	−0.84%

), with an average relative error of only 1.90%, indicating an overall accuracy exceeding 95%. Grassland, which functions as an area adjuster in the model, displayed larger deviations in certain years, such as an 8.38% deviation in 2010. Overall, the model exhibits a mean relative error of only 3.79% across various land use categories, indicating a high level of precision.

3.1.2. Prediction of Quantity Structure of Land Use Under Different Scenarios

The SD model sets 2020 as the base year, with a time step of 5 years. After running the model, the forecasted land use quantity structure for Hubei Province from 2020 to 2060 was obtained under three different scenarios (see

Table 4. Land use quantity structure in different SSP-RCP model key years in Hubei Province (unit: km²).

Year	Climate Scenarios	Cultivate Land	Woodland	Grassland	Water	Construction Land	Unused Land
2030	SSP1-1.9	64291.40	91,742.70	5534.16	12,836.00	11,256.81	348.06
	SSP2-4.5	64,328.60	92,238.00	4853.00	12,861.40	11,380.07	348.06
	SSP5-8.5	64,093.30	92,590.40	4784.16	12,700.70	11,492.51	348.06
2060	SSP1-1.9	60,600.50	86,012.10	9995.58	12,540.90	16,475.97	384.08
	SSP2-4.5	60,597.20	89,727.50	5907.55	12,538.60	16,854.20	384.08
	SSP5-8.5	59,794.70	84,273.90	124,12.80	11,990.50	17,153.15	384.08

).

3.2. Simulation of Land Use Structure Under Different Scenarios

3.2.1. Performance Evaluation of the PLUS Model

To validate the training and simulation accuracy of the model, LUCC distributions for 2018 and 2020 were simulated and compared with corresponding observational data. The results demonstrated consistent trends in land use changes, with minor discrepancies observed primarily in peri-urban areas surrounding major cities such as Wuhan and Yichang (see Figure 6).

The validation module in PLUS was used to compare simulated data with observational data for both 2020 and 2018. The kappa coefficients for 2020 and 2018 were 0.9185 and 0.8053, respectively, with an average of 0.8619 over the two periods. The Figure of Merit (FoM) values were 0.1251 for 2020, 0.1103 for 2018, and an average of 0.1177. These metrics indicate a high level of model accuracy.

3.2.2. LUCC Simulation Under Different Scenarios

In this study, conversion rules were applied according to different SSP-RCP pathways. After running the model, LUCC simulation results were obtained for the years 2030, 2040, 2050, and 2060, across three pathways for each period, resulting in a total of 12 maps (Figure 7). These results will be incorporated as factors into stochastic forest models to simulate future plant carbon sinks.

3.3. Simulation of NPP Under Different Scenarios

3.3.1. Evaluation of Random Forest Model

During hyperparameter optimization, the vegetation NPP for 2010 was used as the dependent variable, while the independent variables included vegetation NPP from 2000, land use and average temperature for 2010, and average annual precipitation, as well as static data such as slope, longitude, latitude, and elevation. The RandomizedSearchCV function from scikit-learn was employed for parameter optimization, with the following parameter ranges: n_estimators as random integers from 0 to 800, max_depth selected from [None, 10, 20, 30], min_samples_split ranging from 2 to 10, min_samples_leaf ranging from 1 to 10, and max_features as random integers from 1 to 8. The model was iterated 100 times, with performance evaluated based on mean squared error and R². A 10-fold cross-validation was used for validation, and the results showed that the model performed optimally with n_estimators = 300, min_samples_split = 2, min_samples_leaf = 1, max_features = 8, and max_depth = 50.

From

Table 5. Ten-fold cross-validation results of random forest model.

Folds	Mean Squared Error	R2	Folds	Mean Squared Error	R2
1	1454.76	0.9570	7	1406.97	0.9579
2	1425.71	0.9572	8	1446.69	0.9563
3	1438.60	0.9567	9	1423.14	0.9570
4	1399.83	0.9576	10	1472.63	0.9554
5	1436.41	0.9567	Mean	1435.10	0.9568
6	1446.33	0.9568

, it is evident that the mean squared error (MSE) ranges between 1399.83 and 1472.63, indicating relative stability among the models with only minor random errors. The average R² value is notably high at 0.9568, suggesting a strong model fit, making these models well-suited for subsequent simulation and prediction tasks.

Comparisons across the 30 experiments revealed that while the models performed better during self-cross-validation for 2010, with R² values ranging from 0.958 to 0.960, their performance slightly declined when simulating and observing 2020 data, where R² values ranged between 0.85 and 0.86. The overall relative error increased from 0.89% to 0.95% in 2010 to between 1.55% and 2.08% in 2020. Nevertheless, the average R² remained high at 0.854 during the 2020 validation, with an average relative error of only 1.73%. Scatter and trend analysis (Figure 8) revealed a fitting slope of 0.98 ± 0.01 and a Pearson’s r correlation coefficient of 0.92, indicating high overall model accuracy.

While the models displayed some variability in performance metrics, they generally demonstrated stability, confirming their suitability for subsequent research endeavors.

Figure 9 compares the spatial distribution of MOD17A3 NPP data with simulated values for Hubei Province in 2020. The overall spatial distribution of the simulated values closely aligns with that of the MOD17A3 NPP data. The study also focused on specific regions, such as the eastern plains and western mountainous areas, where the spatial distributions of MOD17A3 NPP and the simulated values were found to be consistent across these distinct geographical settings. This comparative analysis highlights the model’s high fitting accuracy and minimal errors, demonstrating its strong applicability for simulating vegetation NPP in Hubei Province.

3.3.2. NPP Simulation Under Different Scenarios

The model simulations provided NPP projections for Hubei Province for the years 2030, 2040, 2050, and 2060 under three different scenarios. The results indicate that the SSP1-1.9 scenario yields the highest carbon sequestration, while the SSP5-8.5 scenario shows the lowest. The annual carbon sequestration totals range from 115.99 TgC to 117.59 TgC, averaging 116.89 TgC across the three scenarios, representing a growth of 10.16% to 11.68% compared with 2020 levels.

As illustrated in Figure 10, the spatial patterns of future carbon sequestration in Hubei Province are generally consistent across the three scenarios, characterized by a gradient of increasing values from east to west, with higher carbon sequestration levels in the western regions and lower levels in the eastern regions.

4. Discussion

4.1. Feasibility of Predicting Future Spatial Distribution of Carbon Sinks Using Correlation Models

Theoretically and practically, it is feasible to predict the future spatial distribution of carbon sinks based on correlation models. To explore the spatial patterns of future carbon balance, we developed the SD-PLUS model to predict land use quantity structures and simulate spatial distributions under different SSP-RCP scenarios. Additionally, a random forest model was constructed to simulate the spatial distribution of future carbon sequestration. Model accuracy assessments revealed that the SD model has an average relative error of only 3.79% across various land use categories. The PLUS model achieved average kappa and FoM coefficients of 0.8619 and 0.1177, respectively. The random forest model demonstrated high accuracy in both training and simulation, with average R² values of 0.9585 and 0.8541 and relative errors within 2.08%. The simulated NPP values show strong spatial alignment with observed distributions. These results highlight the excellent performance and high precision of the models during both the training and testing phases, affirming their effectiveness in analyzing carbon balance patterns in Hubei Province under different SSP-RCP pathways. This confirms the feasibility of this approach for research purposes.

4.2. Regional LUCC Changes Under Different SSP-RCP Scenarios

In terms of spatial distribution, the future functional zoning of Hubei Province converges across the three scenarios. The western region is predominantly focused on ecological land use, while the central region integrates agricultural production with urban development. The primary land use types exhibit the following characteristics:

The area of construction land increases in all three scenarios, primarily concentrated in the eastern plains of Hubei, particularly within the Wuhan metropolitan area. Notably, SSP5-8.5 shows greater interconnectivity among cities within the Wuhan metropolitan area than SSP1-1.9. In the western regions, such as Enshi Prefecture and Shiyan City, there is minor peripheral expansion. Figure 11 indicates that construction land primarily originates from the conversion of arable land and forests.

In the eastern region, arable land is mainly converted, especially on the outskirts of urban areas within the Wuhan metropolitan area, indicating significant pressure on arable land conservation in these regions. However, due to existing Chinese policies on arable land protection, which require compensatory conversion from non-arable land to arable land when arable land is converted to other types, the loss of arable land in the eastern region is compensated for mainly in the western region. In these areas, arable land is primarily converted to forests and construction land. This results in significantly lower arable land loss under SSP1-1.9 compared with SSP5-8.5.

Forests primarily experience conversion, with significant loss concentrated in the western region. The main conversions include replenishing arable land and transitioning to construction land. The western region serves as the main carbon sink area in Hubei Province, and any reduction in forest land would adversely affect carbon neutrality efforts. SSP5-8.5 shows notably higher forest land loss compared to the other two scenarios. Water bodies, grasslands, and unused land constitute a smaller proportion of LUCC in Hubei Province. Water bodies have a relatively fixed spatial distribution and a limited impact on future carbon balance patterns.

4.3. Variation Trends of NPP Under Different SSP-RCP Scenarios

In Hubei Province, the total NPP follows three distinct curve patterns under the different pathways: “inverse S-shaped”, “S-shaped”, and “parabolic” (Figure 12). Under the SSP1-1.9 pathway, the NPP trend exhibits an “inverse S-shaped” curve, characterized by rapid growth in the early phase, a slower increase in the middle phase—likely due to lower temperatures and reduced rainfall—and a renewed rapid growth in the later phase. This pattern can be attributed to initial rapid growth driven by existing policies, a subsequent slowdown influenced by climatic factors, and a final surge due to the preservation of significant ecological areas, such as forests and grasslands.

The SSP2-4.5 pathway exhibits a “parabolic” trend. This may be because SSP2-4.5 closely reflects current development trends, resulting in a situation where, after prolonged development, vegetation in the region stabilizes and later experiences insufficient growth momentum, reaching a peak. The “parabolic” trend observed in the NPP changes over the long-term sequence reflects this pattern.

The SSP5-8.5 pathway exhibits an “S-shaped” trend in NPP changes. During the initial and middle periods, NPP experiences rapid growth due to increased temperatures that enhance plant photosynthetic efficiency and the CO₂ fertilization effect from rising CO₂ concentrations. However, in the later stages, agricultural expansion encroaches on forest land, leading to deforestation in the western regions. Additionally, temperature increases negatively impact NPP in parts of eastern Hubei. These factors combined to result in sluggish carbon sequestration in the later stages of SSP5-8.5.

Spatially, the Jianghan Plain and the southeast riverine plains of Hubei Province consistently exhibit lower vegetation NPP compared to other regions. These areas are predominantly agricultural, cultivating seasonal crops such as rice and wheat, which have unstable ecological stability and are highly influenced by human activities. By contrast, the western and southern regions of Hubei are mountainous, with dense vegetation cover preserving subtropical evergreen and deciduous broad-leaved mixed forests. Higher temperatures in these areas extended the plant photosynthesis periods, enhancing carbon accumulation.

When comparing the carbon sink distribution in 2020 with projections for 2060, high NPP regions in Hubei Province have gradually shifted northward. Under SSP1-1.9 scenarios, by 2060, carbon sink areas have moved from northwestern Hubei to southwestern Hubei. For instance, in Enshi Prefecture, the average NPP in 2060 remains comparable to that in 2020 at 716.22 gC/(m² × a), whereas Shiyan City in the northwest shows significantly increased NPP in 2060 at 829.12 gC/(m² × a), compared with 704.37 gC/(m² × a) in 2020. This shift is attributed to improved thermal conditions favoring plant growth in the higher latitude regions of northwestern Hubei, thereby accelerating forest growth and enhancing carbon sequestration capacity. Such changes could lead to a negative impact on traditional agriculture in the region, resulting in a decline in agricultural output and a reduction in related employment opportunities. Additionally, climate change may threaten biodiversity, disrupting the original ecological balance. In the southern regions, issues such as soil erosion and land degradation could be further exacerbated. To address these challenges, Hubei Province should strengthen the implementation of ecological protection policies, promote the development of a green economy, raise public awareness of ecological conservation, and guide land use and industrial layout through scientific planning.

4.4. Developing Scientifically Sound Land Use Management Policies Is Crucial for Achieving Effective Carbon Neutrality

Under the SSP1-1.9 pathway, Hubei Province is better positioned to achieve carbon neutrality, with a future capacity to support approximately 4.8 million tons of standard coal in fossil fuels. By controlling the carbon emission coefficient of construction land, effective monitoring and management of land carbon emissions can be ensured, providing a scientific basis for formulating and adjusting land management policies. For instance, maintaining the carbon emission coefficient of construction land between 6800 to 7100 tons/km² per year will likely facilitate carbon neutrality in Hubei Province. Integrating carbon sequestration functions into the “dual evaluation” process for optimizing national land space and implementing “one-loss-one-compensation” and “one-loss-multi-compensation” land policies will effectively ensure that total carbon sequestration remains unchanged during land use changes. By reasonably delineating carbon sequestration functional zones based on future carbon sink distribution—such as forestry carbon sequestration protection zones, forestry carbon sequestration enhancement zones, farmland carbon sequestration enhancement zones, and urban carbon sequestration enhancement zones—targeted “enhancement” strategies can be proposed. These strategies can help implement cross-regional carbon sequestration compensation measures, further exploring the potential of carbon sinks and alleviating carbon neutrality pressures.

5. Conclusions

This study demonstrates that machine learning algorithms, such as the SD-PLUS model and random forest model, can effectively incorporate climate change, LUCC, and socioeconomic development scenarios into the prediction and simulation of NPP. The simulations exhibit high accuracy, indicating that the results can provide reliable data support for government carbon neutrality policies. The application of these machine learning methods in a region with a complex ecosystem in central China reveals their applicability and reliability in simulating NPP trends under different scenarios. The SD-PLUS model describes LUCC as an organic, spontaneous patch generation process, capable of simulating different scenarios of LUCC by coupling various socioeconomic development scenarios. The random forest model successfully constructs relationships between driving factors and NPP, and optimizing critical model parameters through random exploration effectively captures historical NPP variations for future predictions.

Moreover, this study demonstrates that changes in carbon sequestration rates vary under different pathways. For instance, NPP changes in the study area exhibit “inverse S-shaped”, “S-shaped”, and “parabolic” curves under different policies, indicating significant impacts on local ecosystems, and the main carbon sink areas in the region show a northward shift. In scenarios with low carbon emissions, future NPP is projected to have greater growth potential, making it easier to achieve carbon neutrality goals.

However, this study only considers vegetation carbon sinks, neglecting other pathways such as geological and artificial carbon sinks. Additionally, factors influencing vegetation carbon sinks, such as soil properties and soil erosion, were not discussed in the NPP modeling process. These aspects will be addressed in future research endeavors.