Assessing Spatiotemporal Dynamics of Net Primary Productivity in Shandong Province, China (2001–2020) Using the CASA Model and Google Earth Engine: Trends, Patterns, and Driving Factors

Abstract

Net primary productivity (NPP) is a core ecological indicator within terrestrial ecosystems, representing the potential of vegetation growth to offset anthropogenic carbon emissions. Thus, assessing NPP in a given region is crucial for promoting regional ecological restoration and sustainable development. This study utilized the CASA model and GEE to calculate the annual average NPP in Shandong Province (2001–2020). Through trend analysis, Moran’s Index, and PLS−SEM, the spatiotemporal evolution and driving factors of NPP were explored. The results show that: (1) From 2001 to 2020, NPP in Shandong showed an overall increasing trend, rising from 254.96 to 322.49 g C·m⁻²/year. This shift was accompanied by a gradual eastward movement of the NPP centroid, indicating significant spatial changes in vegetation productivity. (2) Regionally, 47.9% of Shandong experienced significant NPP improvement, 27.6% saw slight improvement, and 20.1% exhibited slight degradation, highlighting notable spatial heterogeneity. (3) Driver analysis showed that climatic factors positively influenced NPP across all four periods (2005, 2010, 2015, 2020), with the strongest impact in 2015 (coefficient = 0.643). Topographic factors such as elevation and slope also had positive effects, peaking at 0.304 in 2015. In contrast, human activities, especially GDP and nighttime light intensity, negatively impacted NPP, with the strongest negative effect in 2010 (coefficient = −0.567). These findings provide valuable scientific evidence for ecosystem management in Shandong Province and offer key insights for ecological restoration and sustainable development strategies at the national level.

1. Introduction

Nature is universally acknowledged as the bedrock of human existence and progress [1], bestowing upon us a tapestry of diversity and intricacy that furnishes humanity with a cornucopia of resources and a habitat for life. Vegetation, a pivotal element within the terrestrial ecosystems of our natural world, serves not only as the provider of essential habitats and sustenance for our survival but also stands as a sentinel in the battle against climate change and a key regulator in the global carbon cycle [2]. NPP, as a measure of net carbon sequestration, can be utilized to evaluate the potential of vegetation growth in offsetting anthropogenic carbon emissions [3]. The spatiotemporal distribution of NPP affects energy flow and nutrient cycling within ecosystems, which in turn impacts biodiversity and the stability of ecosystems. A decrease in NPP implies a reduction in the ecosystem’s carbon absorption capacity, exacerbating the greenhouse effect, and could lead to species extinction and ecological degradation. Furthermore, changes in NPP have a direct impact on agricultural productivity, influencing global food production and food security [4]. Additionally, natural factors (such as population distribution, droughts, and floods) [5], and human activities (urbanization, agricultural development, and land management policies) interact to jointly affect the productivity levels of ecosystems [6,7]. An in-depth investigation into the mechanisms driving NPP changes can provide critical insights into how we might address climate change, protect ecosystems, and promote sustainable development.

As public awareness of ecological conservation continues to rise and the concept of sustainable development gains increasing attention, both domestic and international scholars have conducted studies on vegetation NPP and its responses to environmental changes at various scales [8]. As early as the 18th century, researchers began studying NPP to assess forest productivity. Over time, NPP research expanded to aquatic and terrestrial ecosystems [9], where it has been used to quantify the NPP of vegetation. NPP measurement methods can be broadly categorized into two types: field measurements (in situ observations) and model estimations (mathematical simulations) [10]. In the early stages of NPP research (1960s–1970s), studies primarily relied on field measurements. While this traditional method, based on manual, site-specific, and time-fixed measurements, provided accurate and reliable local data, it also had limitations, such as difficulty in extrapolating results to larger regions and challenges in applying these measurements for NPP estimations across different areas. In contrast, model estimations offer a solution to these challenges, enabling effective NPP estimations across broad spatial scales. NPP estimation models can be categorized into three main types: climate productivity models, biogeochemical models, and light use efficiency models. Climate productivity models develop regression relationships between climate factors (such as annual average temperature and precipitation) and field-measured NPP. Common models in this category include the Miami model [11], the Thornthwaite Memorial model [12], and the Chikugo model [13]. While these models are simple and convenient, they may produce significant errors when applied to different study areas [14]. Biogeochemical models estimate NPP by simulating physiological processes, including photosynthesis, organic matter decomposition, and nutrient cycling in plants. Notable models in this category include the BEPS model, the Century model, and the Biome-BGC model. With the advancement of remote sensing (RS) technology, the 1990s saw the application of RS models for NPP estimation. These models utilize satellite imagery to capture large-scale surface data, driving a shift toward NPP estimation over broad spatial scales. Light use efficiency models estimate NPP by utilizing RS data and calculating it based on plant photosynthesis and light use efficiency. Prominent models in this category include the Carnegie–Ames–Stanford approach (CASA) model and the GLO-PEM model [15]. The CASA model estimates NPP using the normalized difference vegetation index (NDVI) and incorporates climate data, including temperature, precipitation, and solar radiation. Its advantages include requiring a relatively small number of RS input parameters, which minimizes potential errors and allows for dynamic monitoring of NPP across various regions and time periods. The model has played a significant role in vegetation productivity simulations across many regions worldwide [16,17]. Large-scale, long-term model estimations pose new challenges to computational power, but the advent of cloud computing has significantly reduced the time costs of NPP estimation over extended periods and large areas. Google Earth Engine (GEE), a powerful cloud computing platform, allows researchers to rapidly process and analyze RS data on a global scale, greatly improving the efficiency of NPP estimation. In recent years, it has become a mainstream method for NPP studies. For example, Ji et al. [18] utilized GEE in combination with the CASA model to monitor and accurately assess NPP in the grasslands of Inner Mongolia over a 20-year period. Furthermore, the advantages of GEE are also evident in its ability to more effectively estimate NPP in highly urbanized regions and complex terrains. Wu et al. [19] improved the CASA model by integrating multi-source RS data, enhancing the accuracy of NPP estimation, particularly in the alpine grasslands of the Qinghai Lake Basin. In summary, the application of RS technology has significantly advanced NPP estimation methods, and the use of the GEE cloud platform has made long-term, large-scale spatiotemporal NPP research feasible. The combination of GEE and the CASA model not only improves computational efficiency but also holds broad application potential on a global scale, making it an ideal tool for large-scale NPP estimation, monitoring, and assessment.

Spatiotemporal evolution analysis can identify and quantify the changing trends of terrestrial vegetation’s physical and chemical parameters. In particular, long-term evolution analysis based on RS satellite data has become a focal point in the field of ecological RS research [20,21,22]. This approach provides continuous, large-scale observational data, offering scientific evidence for understanding ecosystem dynamics and guiding ecological conservation and land use management. For example, Mehmood et al. [23] investigated the spatiotemporal evolution of the NDVI in Khyber Pakhtunkhwa (KPK) Province, Pakistan, from 2000 to 2022, revealing a continuous increase in vegetation cover and its complex adaptation to climate change. L. Liu et al. [24] studied the spatiotemporal evolution of grassland fraction of vegetation coverage (FVC) from 1982 to 2021, uncovering the seasonal impacts and lag effects of various climatic factors on grassland coverage. Hou et al. [25] analyzed the spatiotemporal evolution of NPP in the UANSTM urban cluster on the northern slope of the Tianshan Mountains from 2001 to 2020, quantifying the contributions of aridification and climate change to NPP variations. These studies emphasize the critical role of spatiotemporal analysis in exploring vegetation dynamics, evaluating ecological restoration efforts, and informing environmental management decisions. With long-term RS data, scientists can more accurately monitor and predict changes in vegetation, drought, productivity, and ecosystem services, while assessing the impacts of various ecological projects. This provides crucial scientific evidence for global climate change research and sustainable development. As a key indicator of vegetation growth and ecosystem productivity, NPP is influenced by multiple factors, including climate change, land use/land cover (LULC), and human activities. Previous studies have highlighted the profound impact of LULC on global ecosystems and primary productivity [26,27], offering valuable insights into the ecological consequences of these transformations. Moreover, scientists have extensively explored how human activities, such as fossil fuel combustion and LULC changes, influence the global carbon cycle and the dynamics of primary productivity [28]. Scholars have also assessed the performance of various ecological models in simulating NPP under the context of changing age structures, further revealing the interactions between climate change, LULC, and ecosystem functions. Although extensive research has been conducted on the relationship between NPP and its driving factors, much of this research has relied on linear correlation tests or regression methods, which struggle to capture the complex nonlinear relationships among natural factors. In recent years, many scholars have investigated vegetation NPP and its driving mechanisms [29,30]. Given that vegetation NPP does not respond to different driving factors in a simple linear “driver-response” manner, but rather involves complex multi-scale (both temporal and spatial) interaction processes and response characteristics, it is crucial to establish a nonlinear analytical framework. Adopting models such as structural equation modeling (SEM) has become particularly important for comprehensively evaluating NPP’s nonlinear responses to climate change, LULC, and human activities [31]. Compared to traditional linear models, SEM can account for interactions among multiple variables and alleviate multicollinearity issues by constructing latent variables, thus more accurately revealing complex causal relationships.

As one of China’s key agricultural production areas, Shandong Province has long been regarded as a major growing region for important crops such as wheat, corn, and peanuts. The region experiences a humid climate and sparse surface vegetation, making it one of the areas most significantly affected by global climate change. Over the past 30 years, the average temperature in Shandong has risen by 1.2 °C due to global warming, while annual precipitation has exhibited an irregular fluctuation pattern, leading to variations in vegetation productivity. In addition, long-term human activities, such as urbanization, deforestation, and land use policies, have resulted in significant LULC in Shandong Province. Data [32] show that over the past 20 years, impervious surfaces in the province increased by 48.5%, forest cover expanded by 15.8%, while cultivated land decreased by 9.4%. These changes directly impact the region’s carbon balance and the NPP of vegetation. Recent reports have also indicated that climate change and human activities may lead to vegetation degradation in Shandong Province, potentially triggering the risk of secondary desertification. However, most of these studies are limited to localized observations or surveys, and lack in-depth analysis of the overall vegetation trends and driving factors for the province. Moreover, previous research has largely focused on either climate change or LULC when assessing their impacts on vegetation NPP, failing to comprehensively evaluate the relative contributions of climate factors, natural environmental factors, and human activities to NPP variation. Therefore, it is crucial to investigate the spatiotemporal changes of vegetation NPP and their driving mechanisms in Shandong Province through modeling studies.

Building on this background, a comprehensive study of the mechanisms and influencing factors of NPP spatiotemporal dynamics can offer critical insights into addressing climate change, protecting ecosystems, and promoting sustainable development. This study aims to estimate the NPP of vegetation, reveal its spatiotemporal evolution characteristics, and analyze the various factors influencing its changes. Taking Shandong Province as a case study, the main objectives are as follows: (1) to calculate NPP using the CASA model implemented on the GEE platform; (2) to examine the spatiotemporal dynamics and clustering patterns of NPP in Shandong Province using trend analysis and Moran’s I; and (3) to investigate the comprehensive impacts and relative contributions of 11 driving factors, grouped into climate change, natural factors, and human activities, on NPP changes. This research provides a new perspective for understanding the regulatory mechanisms of regional vegetation productivity. The findings enhance the understanding of the spatiotemporal dynamics of vegetation productivity in Shandong Province and provide a vital theoretical foundation and practical guidance for ecological improvement and sustainable development, both within the region and globally.

2. Methodology and Materials

2.1. Research Area

Farmland dominates approximately 65% of the landscape of Shandong Province (34°23′N–38°24′N, 114°47′E–122°43′E). The region features rainy, humid summers and cold, dry winters. Interannual fluctuations in NPP are primarily driven by seasonal variations in precipitation. On the other hand, with the rapid urbanization of Shandong Province, NPP in certain areas has significantly declined. Additionally, the interaction between changes in the natural environment and socioeconomic development further affects NPP performance across different regions of Shandong Province [33,34]. As a major agricultural province, the NPP of vegetation in Shandong is not only directly related to the regional ecological balance but also indirectly linked to carbon sequestration and regional climate regulation functions. Therefore, under the context of global climate change and ecological protection, studying NPP in Shandong Province is of great significance for ecological conservation, climate change mitigation, resource management, and socioeconomic development. On the one hand, NPP is a crucial indicator of vegetation growth and ecosystem productivity, playing a key role in assessing regional carbon balance, climate change, and ecological health. On the other hand, a comprehensive analysis of the spatiotemporal variations in NPP and their driving factors can offer valuable scientific insights for developing LULC policies and ecological protection strategies. Given the regional differences in natural environments and urbanization intensity, this study divides Shandong Province into five subregions (Figure 1).

2.2. Data Source

The datasets are summarized in

Table 1. Overview of data sources, spatial and temporal resolutions, and periods used in the study.

Data	Abbreviation	Spatial Resolution	Temporal Resolution	Period	Data Source
Normalized Difference Vegetation Index	NDVI	250 m	16-day	2001–2020	GEE: MODIS/006/MOD13Q1
ECMWF ERA5-Land Solar Radiation	SOL	9 km	Daily		GEE: ECMWF/ERA5_LAND/DAILY_AGGR
MODIS Evapotranspiration	ET, PET	500 m	8-day		GEE: MODIS/006/MOD16A2
Average Annual Precipitation	PRE	1 km	Monthly		[35]
Annual Average Temperature	TEM		Monthly		[35]
NASA SRTM Digital Elevation	DEM	30 m	N/A	2000	GEE: USGS/SRTMGL1_003
Soil Type	Soil	1 km	Annual	2023	FAO: Harmonized World Soil Database
Gross Domestic Product	GDP			2001–2020	[36]
China Land Cover Dataset	LAND				[32]
Nighttime Light (2001–2013)	Light			2001–2013	GEE: DMSP/OLS Nighttime Lights (NOAA/DMSP-OLS/NIGHTTIME_LIGHTS)
Nighttime Light (2014–2020)	Light	500 m		2014–2020	GEE: VIIRS Nighttime Day/Night Band
Population Density	POP	100 m		2001–2020	GEE: WorldPop (WorldPop/GP/CHI_POP)
MODIS Land Surface Temperature	LST	1 km	8-day		GEE: MODIS/006/MOD11A2

. The input variables for the CASA model comprise the NDVI, land surface temperature (LST), solar radiation (SOL), evapotranspiration (ET), and potential evapotranspiration (PET). Based on previous studies [7] and the availability of data in the study area, the driving factor data are classified into three main categories: climate factors, natural environment data, and human activity data. Climate factors include annual PRE, annual mean TEM, and PET, which are employed to assess the effects of climate change on NPP. Natural environment data consist of the DEM from NASA, employed to investigate the influence of terrain features such as slope, aspect, and elevation on NPP. Soil type, obtained from the HWSD, is used to evaluate the impact of various soil types on vegetation productivity. Human activity data encompass nighttime light data; GDP data, utilized to explore the relationship between economic development and vegetation productivity; LULC type (we referred to the relevant literature to classify land use types into six LULC categories based on the China land cover dataset: cropland, forest, shrubland, grassland, water bodies, barren land, and impervious surfaces); and population density data, used to assess the population distribution’s impact on NPP. All datasets were spatially processed and reprojected into a unified coordinate system to ensure the robustness and scientific accuracy of the multi-source data fusion analysis. Referring to the relevant literature [7,25], we conducted spatiotemporal analyses and driving factor assessments of the study area at five-year intervals for the years 2001, 2005, 2010, 2015, and 2020.

3. Research Framework and Methods

3.1. Research Framework

The research framework employed in this study is shown in Figure 2. First, the annual average NPP was calculated using the CASA model implemented within the GEE platform. Sen’s slope estimation and the Mann–Kendall trend test were applied to analyze the long-term NPP trends. Second, Global Moran’s I statistics were utilized to assess the spatial autocorrelation of NPP, revealing spatial clustering patterns. Finally, PLS−SEM was applied to identify the key factors driving changes in NPP, including climate factors (temperature, precipitation, evapotranspiration), human activities (GDP, LULC type, nighttime light intensity, population density), and natural environmental factors (slope, aspect, elevation, soil type). This study not only deepens the understanding of vegetation growth dynamics in Shandong Province but also offers valuable theoretical and practical insights for ecological enhancement and sustainable development at both regional and global levels.

3.2. Methods

3.2.1. Estimation of NPP

NPP represents the amount of carbon sequestered by plants from the atmosphere through photosynthesis, serving as a crucial indicator of ecosystem productivity. NPP reflects the efficiency of plant growth and its response to environmental conditions, making it a key parameter for assessing ecosystem health and carbon cycling processes. In this study, long-term NPP in Shandong Province was achieved by the synergistic combination of the CASA model and GEE [37]. The CASA model computes NPP by multiplying the absorbed photosynthetically active radiation (APAR) with the light use efficiency (ε), using the following formula:

$$NP{P}_{(x,t)}=APA{R}_{(x,t)}\times {\epsilon }_{(x,t)}$$

(1)

where $NP{P}_{(x,t)}$ represents the NPP of pixel x at a specific time t, $APA{R}_{(x,t)}$ denotes the solar radiation absorbed by the canopy, and ${\epsilon }_{(x,t)}$ is the actual light use efficiency.

Additionally, APAR is calculated from the NDVI and total solar radiation data, while ε is adjusted based on environmental factors such as temperature and precipitation. The formula is as follows:

$$APA{R}_{(x,t)}={PAR}_{(x,t)}\times FPA{R}_{(x,t)}\times 0.5$$

(2)

where ${PAR}_{(x,t)}$ is the total solar radiation for pixel x at a specific time t, extracted from the short-wave radiation in the ERA5 reanalysis data, and $FPA{R}_{(x,t)}$ represents fraction of photosynthetically active radiation and can be estimated by extracting the NDVI from MOD17A3 products.

The actual light use efficiency (ε) is determined using the following equation:

$${\epsilon }_{(x,t)}=T_{\epsilon 1(x,t)}\times T_{\epsilon 2(x,t)}\times W_{\epsilon (x,t)}\times {\epsilon }_{max}$$

(3)

where $T_{\epsilon 1(x,t)}$ and $T_{\epsilon 2(x,t)}$ represent the stress induced by the temperature of pixel x and specific time t. The detailed algorithms for $T_{\epsilon 1(x,t)}$ and $T_{\epsilon 2(x,t)}$ can be found in [10]. $W_{\epsilon (x,t)}$ represents moisture stress, and ${\epsilon }_{max}$ represents the maximum light use efficiency of a vegetation type and was adopted the approach of Zhu et al. [38].

$$W_{\epsilon (x,t)}=0.5+0.5\times {\displaystyle \frac{{ET}_{(x,t)}}{{PET}_{(x,t)}}}$$

(4)

where $ET$ and $PET$ represent the estimated and potential evapotranspiration (mm), respectively. In this model, $W_{\epsilon (x,t)}$ generally ranges from 0.5 to 1 for very arid ecosystems and very wet ecosystems, respectively [39].

3.2.2. Sen’s Slope Estimation

Sen’s slope is a non-parametric statistical method used to estimate the rate of change or trend in time series data [23]. It calculates the slope between every pair of points in the time series and uses the median of these slopes to determine the overall trend. Unlike traditional linear regression, Sen’s slope does not rely on assumptions about data distribution, making it especially suitable for analyzing time series with irregularities or non-normal distributions. This method effectively quantifies the changing trends of variables over time and is highly robust against outliers. In this study, Sen’s slope was applied to evaluate the interannual variation trends of long-term NPP in Shandong Province. The formula is as follows:

$$Slope=Median\left[{\displaystyle \frac{(X_j-X_i)}{j-i}}\right],\forall 1\le i\le j\le n$$

(5)

where $X_i$ and $X_j$ are the pixels value of the variables at time $i$ and $j$, and $n$ is the length of the time series. The slope represents the trend within the time series data. When the slope (β) is greater than 0, an increasing trend is observed in the NPP, indicating an improvement in vegetation productivity. Conversely, when β is less than 0, the NPP exhibits a decreasing trend, reflecting potential vegetation degradation. If β equals 0, it indicates a stable trend in NPP, suggesting little to no change over the analyzed period. Supplementary S1 provides an example.

3.2.3. Mann–Kendall Test

The Mann–Kendall test is a non-parametric statistical method employed to detect monotonic trends, either consistently increasing or decreasing, in time series data [40]. This test is highly effective in detecting significant trends in time series data and is well-suited for non-normal distributions, as it is robust against the presence of outliers. In this study, the Mann–Kendall test was used to analyze the long-term trends of NPP across various regions of Shandong Province. By applying this test, the significance of NPP trend changes can be assessed, thereby ensuring the statistical reliability of the observed trends. This approach not only confirms NPP change patterns but also validates the results of Sen’s slope analysis. The formula for the Mann–Kendall test is detailed in Supplementary S1.

3.2.4. Global Moran’s I

Global Moran’s I is a statistical index used to evaluate the overall spatial autocorrelation of a dataset, helping to identify whether the data exhibit a random distribution, clustering, or dispersion across the study area [41]. The Moran’s I index measures the degree of spatial clustering, aiding in the identification of spatial patterns, such as areas with high or low values of a given variable. In this study, Global Moran’s I was utilized to examine the spatial distribution patterns of NPP in Shandong Province. By calculating the Moran’s I index, the spatial autocorrelation of NPP can be assessed, highlighting trends of spatial clustering or dispersion. This analysis reveals the spatial structure characteristics of NPP across various regions in Shandong Province, offering valuable insights into the spatial heterogeneity of NPP. The formula for Global Moran’s I is as follows:

$${Moran}^{\mathrm{’}}s I=\frac{N{\sum }_i {\sum }_j w_{ij}(x_i-\mu )(x_j-\mu )}{({\sum }_i {\sum }_j w_{ij}){\sum }_i (x_i-\mu {)}^2}$$

(6)

where $x_i$ and $x_j$ represent the NPP values at i and j, $\mu$ is the average NPP value, $w_{ij}$ denotes the spatial weight between i and j, and $N$ is the total sample size.

3.2.5. Center of Gravity Migration Model in NPP

To analyze the spatial evolution of NPP in Shandong Province from 2001 to 2020, this study utilizes the gravity center migration model [42,43]. This model effectively depicts changes in the spatial distribution of NPP over the study period. The model is expressed as follows:

$${\mathrm{X}}_{\mathrm{t}}=\sum _{\mathrm{i}=1}^{\mathrm{n}} ({\mathrm{C}}_{\mathrm{t}\mathrm{i}}\times {\mathrm{X}}_{\mathrm{t}\mathrm{i}})/\sum _{\mathrm{i}=1}^{\mathrm{n}} {\mathrm{C}}_{\mathrm{t}\mathrm{i}}$$

(7)

$${\mathrm{Y}}_{\mathrm{t}}=\sum _{\mathrm{i}=1}^{\mathrm{n}} ({\mathrm{C}}_{\mathrm{t}\mathrm{i}}\times {\mathrm{Y}}_{\mathrm{t}\mathrm{i}})/\sum _{\mathrm{i}=1}^{\mathrm{n}} {\mathrm{C}}_{\mathrm{t}\mathrm{i}}$$

(8)

where ${\mathrm{X}}_{\mathrm{t}}$ and ${\mathrm{Y}}_{\mathrm{t}}$ represent the longitude and latitude coordinates of the NPP gravity center in year $\mathrm{t}$, respectively; $\mathrm{n}$ represents the total number of patches in year $\mathrm{t}$; ${\mathrm{C}}_{\mathrm{t}\mathrm{i}}$ denotes the NPP value of the $\mathrm{i}$-th patch in year $\mathrm{t}$; ${\mathrm{X}}_{\mathrm{t}\mathrm{i}}$ and ${\mathrm{Y}}_{\mathrm{t}\mathrm{i}}$ are the geometric center coordinates of the $\mathrm{i}$-th patch in year $\mathrm{t}$.

3.2.6. Correlation Analysis

In this study, the correlation coefficient method was used to calculate the relationships between NPP and various driving factors, including climatic factors, natural environmental factors, and human activity factors, to analyze their spatial and temporal patterns. The formula for calculating the correlation coefficient between two variables is as follows [44]:

$$r={\displaystyle \frac{\sum _{i=1}^n (x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum _{i=1}^n (x_i-\overline{x}{)}^2\sum _{i=1}^n (y_i-\overline{y}{)}^2}}}$$

(9)

where $n$ is the sample size, $x_i$ represents the observed NPP value at time $i$, $y_i$ represents the observed value of the driving factor at time $i$, and $\overline{x}$ and $\overline{y}$ are the mean values of NPP and the driving factor, respectively. When $r>0$, it indicates a positive correlation between NPP and the driving factor, meaning that as the driving factor increases, NPP also increases. Conversely, when $r<0$, it indicates a negative correlation, meaning that as the driving factor increases, NPP decreases.

3.2.7. The PLS−SEM Model

PLS−SEM is a multivariate analysis method that estimates causal relationships between latent variables while simultaneously constructing and validating these latent variables using observed indicators. PLS−SEM can effectively address multicollinearity issues in complex systems and quantify both the direct and indirect effects of multiple factors on a target variable, making it a crucial tool for analyzing intricate causal relationships [12,45,46]. In this study, PLS−SEM was employed to identify the contributions of climatic conditions, human activities, and natural environmental factors to changes in NPP and to quantify their direct, indirect, and total effects. The conceptualized model of the drivers of NPP in Shandong is presented in Figure 3. The influence of each factor is represented by the path coefficient, with the total effect calculated as the sum of direct and indirect effects. This analysis provides scientific insights into the dynamic changes of Shandong’s ecosystem and helps identify key strategies for ecological management.

To identify the best-fitting model, R², Q², and the goodness of fit (GOF) values were selected as model fit indicators [47]. The goodness-of-fit for each indicator is shown in

Table 2. PLS−SEM evaluation criteria [47].

Criteria	Value	Description
R2	>0.67	Substantial explanatory power
	>0.33	Moderate explanatory power
	>0.19	Weak explanatory power
Q2	>0	A larger value denoting higher prediction accuracy of the model
GOF	0.1	Overall fit of the model is weak
	0.25	Overall fit of the model is medium
	0.36	Overall fit of the model is strong

. The model construction, evaluation, and modification were all performed using the SmartPLS 3 software [48]. Once the model is successfully fitted, the standardized regression coefficients for each path will appear. The direct effect refers to the path coefficient that directly influences NPP changes, while the indirect effect is calculated as the product of path coefficients that impact NPP changes through other variables. The total effect is the sum of both direct and indirect effects.

4. Results

4.1. Spatiotemporal Changes in NPP

4.1.1. Temporal Changes

To uncover the temporal evolution characteristics of NPP in Shandong, linear regression analysis was conducted. The overall interannual (2001–2020) variation in NPP exhibited an upward trend, though the rate of increase is relatively modest in Figure 4. The trend line shows that NPP increased at a rate of 2.0916 g C·m⁻²/year, with an R² value of 0.27, indicating that although the NPP showed an increasing trend over these 20 years, the rate was relatively slow. The annual average NPP remained mostly within the range of 250–350 g C·m⁻²/year during this period, with the lowest value observed in 2002 (224.82 g C·m⁻²/year), and a noticeable peak in 2020 (322.49 g C·m⁻²/year). The overall trend suggests that vegetation productivity in the Shandong region is showing a positive growth trajectory. For different regions within Shandong Province, the NPP trends from 2001 to 2020 exhibited significant variability (Figure 5). In the southwestern region, NPP remained relatively stable throughout the study period, with a slightly upward trend (slope of 0.8755). In contrast, the central−southern and northwestern regions showed more noticeable upward trends, with slopes of 2.071 and 4.5883, respectively, and relatively high R² values, indicating an increase in vegetation productivity over the study period. The northern region also exhibited stable growth, with a slope of 1.2942 and an R² value of 0.1332. In the eastern region, although NPP showed a slight increase (slope of 1.0453), the trend was not significant, with an R² value of only 0.0658. Overall, except for the southwestern region, NPP in other areas of Shandong Province exhibited varying degrees of growth.

Additionally, the shift in the spatial centroid of NPP over the long-term series is represented by centroid migration. The trajectory of NPP centroid migration in Shandong Province from 2001 to 2020 indicates a gradual eastward shift during this period (Figure 6). The NPP centroid migrated from the western region to the eastern region between 2001 and 2020. In 2001, the NPP centroid was located at the position marked by the green triangle, while by 2020, the centroid had shifted to the position marked by the red triangle. This eastward migration of the NPP centroid reflects the spatial distribution changes of NPP throughout the study period, potentially associated with LULC changes and vegetation management practices in the eastern region. Furthermore, the migration path indicates significant spatial changes in vegetation productivity across Shandong Province, with the spatial distribution centroid experiencing a relatively stable, gradual eastward shift over the 20-year period.

4.1.2. Spatial Changes

The spatial distribution of the annual average NPP (2001–2020) in Shandong Province demonstrated significant variation, with higher values predominantly concentrated in the southeast and lower values in the northwest (Figure 7). The trend analysis of NPP (2001–2020) changes revealed pronounced regional differences in spatial distribution and proportions (Figure 7). As outlined in Table 2, NPP trend changes in Shandong Province were categorized into five groups: significant improvement, slight improvement, stable, slight degradation, and severe degradation. Specifically, about 47.9% of the area experienced significant improvement, primarily concentrated in the central and eastern regions of the province, indicating a notable increase in vegetation productivity over the study period. Slight improvement was observed in 27.6% of the area, with a widespread distribution across much of Shandong Province. Areas with stable NPP accounted for just 4.4%, mainly situated in the central-western part of the province, where NPP remained relatively unchanged during the study period. Slight degradation was dispersed across 20.1% of the region, primarily in the southern and western parts of Shandong, indicating a slight decline in vegetation productivity. This suggests that these areas might be under environmental stress or influenced by human activities. Overall, during the 2001 to 2020 period, the spatial distribution of NPP in Shandong Province revealed significant improvement in the central and eastern regions, whereas slight degradation was observed in the southern and western areas. Considering the area proportions of each region, it can be inferred that the overall vegetation productivity in Shandong improved during the study period, though some regions require further attention.

The NPP’s spatial distribution and proportions trend variations in Shandong Province from 2001 to 2020 showed considerable changes across different periods (Figure 8). In 2001–2005, areas with significant improvement in NPP accounted for as much as 91.5%, indicating a substantial increase in vegetation productivity during this period. Areas with slight improvement made up 4.7%, while stable areas accounted for only 3.8%. Areas with slight and severe degradation were negligible, at 0.1% and 0%, respectively. In 2006–2010, the NPP’s spatial distribution changed notably, with the area of significant improvement decreasing to 48.2%, while areas with slight degradation surged to 51.2%. This period reflected the fluctuating nature of vegetation productivity, with signs of degradation emerging in some regions. In 2011–2015, the proportion of areas with significant improvement rebounded to 73.7%, while areas with slight degradation decreased to 22.9%. Stable areas slightly increased to 3.2%, and areas with severe degradation remained minimal at 0.2%. This period indicated a general recovery in NPP across Shandong, though some regions still experienced degradation. In 2016–2020, a significant divergence in NPP trends was observed, with areas of slight degradation increasing sharply to 61.4%, and areas of significant improvement dropping to 38.3%. This period suggests a decline in vegetation productivity in certain regions, potentially due to human activities or environmental pressures. Overall, the NPP trends in the study area from 2001 to 2020 displayed an initial improvement followed by fluctuations. The spatial distribution and proportion changes over the different time periods reflect the complex dynamics of vegetation productivity.

4.1.3. Spatial Autocorrelation of NPP Distribution

The Global Moran’s I value for NPP in 2001, 2005, 2010, 2015, and 2020 were 0.695, 0.723, 0.683, 0.778, and 0.640, respectively. All of these values passed the significance level test, indicating a significant spatial autocorrelation of NPP in these years (Figure 9). The Moran’s I scatter plot is divided into four quadrants: HH (High−High), HL (High−Low), LH (Low−High), and LL (Low−Low), representing four types of local spatial autocorrelation. In 2001, high−value clusters of NPP were mainly located in the eastern and central−southern regions, while low−value clusters were concentrated in the western and northern regions. Spatial autocorrelation in the central and southwestern areas was not significant. By 2005, high−value clusters remained in the eastern and central−southern regions, with low-value clusters persisting in the western and northern parts. In 2010, high−value clusters in the eastern and central−southern regions expanded, while low−value clusters in the western region remained relatively stable. Meanwhile, some areas in the central and southwestern regions exhibited insignificant spatial autocorrelation. In 2015, the central−southern region continued to be a high−value clustering area for NPP, while the southwestern and northern regions were low-value clusters. By 2020, the central-northern region emerged as a high−value cluster for NPP, while the southwestern and northern regions continued to display low−value clusters, with some areas in the central and eastern regions showing insignificant spatial autocorrelation. Overall, the spatial clustering of NPP in Shandong Province exhibited a pronounced clustering effect throughout the study period, particularly in the eastern and central−southern regions, where NPP demonstrated a strong clustering trend. In contrast, the western and northern regions remained low−value clusters, highlighting a spatial disparity in vegetation productivity across the province.

4.2. Driving Factors Analysis for NPP

4.2.1. Correlation Analysis of NPP

The Pearson correlation coefficient was applied for the correlation analysis, and Figure 10 illustrates the relationships between NPP and 13 other variables. Over the four periods of 2005, 2010, 2015, and 2020, the correlation coefficients between different variables exhibited significant temporal variations. In 2005, NPP exhibited a positive correlation with temperature (TEM) and potential evapotranspiration (PET), while it was negatively correlated with precipitation (PRE) and GDP. Elevation (ELE) and slope (SLOPE) also demonstrated a significant negative correlation with NPP. Over time, these relationships evolved. By 2010, the positive correlation between NPP and temperature had weakened, whereas the negative correlation with precipitation and population density (POP) had intensified. Meanwhile, a positive correlation between GDP and other human activity-related variables, such as nighttime light intensity (light), began to emerge, indicating the increasing impact of economic activity on NPP. In 2015, the positive correlation between temperature and NPP strengthened again, and the link between potential evapotranspiration and NPP became more pronounced. However, the negative correlations between NPP, precipitation, and GDP persisted and grew more significant. By 2020, the positive correlations between NPP, temperature, and evapotranspiration remained notable, while the negative correlations with precipitation, GDP, and nighttime light intensity became even more evident. This trend suggests that, over time, the interactions between climate factors and human activities affecting NPP have become increasingly complex. Additionally, the negative correlations between NPP, elevation, and slope remained consistent across all periods, particularly in mountainous and hilly areas, indicating that topographic factors continued to constrain vegetation productivity in the long term. Overall, Figure 8 highlights the complex interactions between NPP and various factors, including climate, natural environmental factors, and human activities, in Shandong Province from 2001 to 2020. These relationships have shown dynamic changes over time, providing valuable insights into the mechanisms driving ecosystem changes in the region.

4.2.2. Structural Equation Modeling (SEM) of NPP Drivers PLS−SEM

An analysis of the factors influencing NPP was conducted using the PLS−SEM method. The model elucidates NPP variation by utilizing path coefficients between different variables, highlighting the direct and indirect effects of climate factors, human activities, and topographic features on NPP. Figure 11 presents the PLS−SEM model results for the distribution of variables and NPP in 2005 (a), 2010 (b), 2015 (c), and 2020 (d). The Q² values for each year were all greater than zero, indicating high predictive accuracy for the model. Additionally, the GOF values also showed that the overall model fit ranged from 0.301 in 2005 to 0.322 in 2020, indicating a moderate level of fit (GOF > 0.25), with the GOF value in 2015 reaching 0.345, approaching a strong level of fit (GOF > 0.36). The structural equation models for the factors influencing NPP in Shandong Province in 2005, 2010, 2015, and 2020 (Figure 11) indicate that changes in human activities had the largest negative impact on NPP, with impact coefficients of −0.499, −0.567, −0.44, and −0.559, respectively. The absolute values of the path coefficients ranged from −1 to 1, with higher absolute values indicating a greater impact of the variable on the target variable. Among the variables, nighttime light and GDP had the highest weights, with values exceeding 0.85 across all four periods, indicating that urban expansion and the degree of urbanization had a considerable negative impact on vegetation NPP. The weight of LULC was approximately 0.5, implying that changes in LULC also significantly affected vegetation NPP. In contrast, population density had a relatively minor impact on NPP.

The total impact of topographic factors on NPP in Shandong Province is shown in

Table 3. Influence coefficients of various factors.

Year	Influencing Factor	Direct Effect	→Climate	→Human	Total Indirect Effect	Total Effect
2005	Climatic Factors	0.349		0.031	0.031	0.38
	Natural Environmental Factors	0.062	−0.032	−0.004	−0.036	0.026
	Human Activity Factors	−0.499	/	/	/	−0.499
2010	Climatic Factors	0.206		−0.009	−0.009	0.197
	Natural Environmental Factors	0.017	−0.009	−0.002	−0.011	0.006
	Human Activity Factors	−0.567	/	/	/	−0.567
2015	Climatic Factors	0.643		0.006	0.006	0.649
	Natural Environmental Factors	0.304	−0.149	−0.038	−0.187	0.117
	Human Activity Factors	−0.44			/	−0.44
2020	Climatic Factors	0.163		−0.0001	−0.0001	0.1629
	Natural Environmental Factors	0.232	−0.105	−0.03	−0.135	0.097
	Human Activity Factors	−0.559			/	−0.559

Note: The arrow (→) in the table indicates the indirect effect of each factor on NPP.

, with values of 0.026, 0.006, 0.117, and 0.097 for the respective years. The main impact pathways include the direct influence of topography on NPP (0.062, 0.017, 0.304, 0.232), the indirect impact through changes in annual climate factors affecting NPP (−0.032, −0.009, −0.149, −0.105), and the indirect impact through influencing human activities (−0.004, −0.002, −0.038, −0.03). Overall, elevation, slope, and aspect all have a positive impact on NPP. Among these, elevation and aspect have the greatest influence on NPP, indicating that the favorable topographic conditions in Shandong provide suitable environments for vegetation growth.

NPP in Shandong Province was positively influenced by climatic factors, with total impact coefficients for the four periods being 0.038, 0.197, 0.649, and 0.1629, respectively (Table 3). The main impact pathways encompass the direct impact of climate factors on NPP (0.062, 0.017, 0.304, 0.232) and the indirect influence through human activities on NPP changes (0.031, −0.009, 0.006, −0.0001). However, when analyzed by period, the impact of climate factors on NPP was particularly strong in 2015, with temperature having a significant impact on NPP in Shandong. As the years progressed, the influence of temperature (TEM) and potential evapotranspiration (PET) on NPP increased, while the impact of precipitation shifted from negative to positive overall. This shift may be linked to a decrease in extreme weather events resulting from climate change and enhanced water resource management.

5. Discussion

5.1. Advantages and Limitations of Using the CASA Model Based on GEE to Calculate NPP

Terrestrial ecosystem productivity is a key indicator for representing the material cycles, energy flows, and ecological functions of terrestrial ecosystems [49]. Analyzing the spatiotemporal distribution and long-term dynamic evolution of ecosystem productivity is a primary approach for comprehensively detecting changes in ecosystem productivity. This is also a crucial step in understanding the mechanisms driving productivity changes. Therefore, NPP, as an important indicator for monitoring ecosystems, is widely used in vegetation productivity monitoring and environmental assessment [50]. In traditional methods for estimating NPP, there are primarily two approaches: direct observation and indirect estimation [51]. Direct observation is known for its high accuracy; however, due to the uneven spatial distribution of observation sites, this method struggles to provide a comprehensive reflection of large-scale NPP and its dynamic changes.

In contrast, indirect estimation methods encompass model-based approaches (including statistical models and parametric models) and RS inversion methods [52]. The model-based approach faces certain limitations in practical applications, such as insufficient accuracy, the need for numerous parameters, and the difficulty of obtaining data, all of which can lead to uncontrollable errors. RS inversion, on the other hand, combines the strengths of RS data, geographic information systems (GIS), and model simulations, and has been widely applied to estimate NPP at both global and regional scales [29,53]. Among these methods, the CASA model, a RS-based light use efficiency model, stands out for its simplicity, requiring fewer parameters and providing easy access to multi-temporal data. These advantages have made it a popular tool for simulating NPP across various ecosystems. The model effectively estimates NPP and describes its dynamic changes at different spatial scales, earning widespread recognition and application by researchers globally. A notable example is the work of F. Yang et al. [54], who utilized the GEE platform and the CASA model to develop a monthly high-resolution NPP dataset for the Qinghai Plateau from 1987 to 2021. This study highlights the significant advantages of the CASA model (particularly when integrated with the GEE platform) in large-scale ecosystem studies, especially for estimating vegetation NPP.

RS image processing is extensively used in various fields, including LULC classification, hydrological studies, urban development strategizing, evaluation of natural disaster impacts, and climatological examinations. Researchers from various disciplines have utilized GEE [55] to tackle their big data processing challenges and continue to make substantial progress in addressing global issues involving geographic big data. With its global data access and powerful cloud computing capabilities, the GEE platform allows researchers to efficiently process and analyze long-term RS data, even without local computing resources. The CASA model, based on light use efficiency and integrating RS data for NPP inversion, provides benefits such as simplicity, fewer parameters, and easy access to multi-temporal data [16,19,56]. It can simulate NPP across various ecosystems and perform NPP inversion and dynamic description at different spatial scales, gaining widespread recognition and application by scholars worldwide. This has been validated in the present study for large-scale NPP estimation in Shandong Province. The combination of the GEE platform with the CASA model has greatly improved data processing efficiency and estimation accuracy through large-scale automated processing and analysis [57]. This is particularly evident in large-scale, long-term ecological studies. However, the RS and meteorological data provided by the GEE platform still have limitations in terms of spatial and temporal resolution. For example, lower-resolution LULC type data can introduce estimation errors and may fail to accurately reflect small-scale surface changes. Additionally, NPP estimates for water bodies, such as rivers and lakes, could encounter inaccuracies stemming from the issues posed by mixed pixels in remotely sensed imagery. To enhance model accuracy, future studies should incorporate higher-resolution data and use ground-based observations for cross-validation. In conclusion, while the CASA model on the GEE platform provides a powerful tool for NPP estimation, users should be mindful of its limitations and make appropriate parameter adjustments and model corrections based on specific conditions.

5.2. Spatiotemporal Evolution of NPP

The study of spatiotemporal dynamics is a core aspect of research on ecosystem productivity and the carbon cycle. By analyzing the spatiotemporal evolution of NPP, insights can be gleaned into the processes through which ecosystems adapt to climate fluctuations and the pressures exerted by human actions. This analysis offers crucial insights for ecological restoration, land management, and climate adaptation at both regional and global scales [58]. Within the contemporary research landscape, the proliferation of RS technology and big data analysis has significantly advanced the study of NPP spatiotemporal changes. At the national level, a multitude of researchers have delved into the intricate interplay linking NPP with LULC dynamics. For example, Fu et al. [59] demonstrated that changes in LULC, such as cropland expansion, urbanization, and forest restoration, significantly influence the spatiotemporal dynamics of NPP. These changes highlight the roles that different land use types play in ecosystems and reflect how human activities influence ecosystem productivity by altering land use patterns.

Shandong, a pivotal epicenter for agriculture and industry within China, has experienced significant transformations in land utilization patterns, which in turn have exerted a substantial influence on the spatiotemporal dynamics of NPP, particularly under the backdrop of rapid urbanization and industrialization. The substitution effect of urban expansion on natural vegetation has significantly influenced NPP’s spatial distribution. NPP in Shandong Province has shown a general increase (254.96 to 322.49 g C·m⁻²/year) as can be seen in Figure 4, at 2.0916 g C·m⁻²/year., Notwithstanding the challenges posed by urbanization and anthropogenic interventions, this suggests that the overall vegetation productivity in Shandong Province has improved (particularly in the central−southern mountainous areas and eastern coastal regions, where the enhancement has been most significant). At the same time, the spatial heterogeneity of NPP is quite evident (Figure 5). High NPP values were predominantly observed across the eastern and southern sectors, with notable concentrations on the Jiaodong Peninsula and the mountainous zones situated in the central−southern parts. These regions, located on the windward slopes of the southeastern monsoon and influenced by the monsoon climate, benefit from favorable water and heat conditions that promote robust vegetation growth. In contrast, lower NPP values are observed in Shandong’s northern and western regions, particularly in areas heavily impacted by urban expansion, where signs of NPP degradation are evident. A substantial portion—nearly 47.9% of the regions—witnessed a marked enhancement in NPP, while 27.6% showed slight improvement and 20.1% experienced slight degradation (Figure 7). These findings indicate that variations in land use patterns across different regions of Shandong Province have directly impacted the spatial distribution of NPP. To thoroughly analyze the primary drivers of land use type change, we utilized relevant data from the year 2020 and conducted a correlation analysis between land use types and external variables, including climate, topography, and the impact of human activities. The correlation figures can be found in Supplementary S1. Furthermore, the analysis of centroid migration (Figure 6) reveals that from 2001 to 2020, the NPP centroid gradually shifted eastward. This migration reflects an increase in vegetation productivity in the eastern regions, particularly around the Jiaodong Peninsula. It also confirms the significant role of land use management policies, such as ecological restoration projects. On a national scale, research shows that LULC changes affect NPP differently in various regions. Ref. [60] pointed out that cropland expansion and forest restoration across the country contributed to increased NPP, while urbanization and infrastructure development had negative impacts. The results from Shandong Province align with this trend. Figure 12 and

Table 4. Mean annual NPP (g·C·m⁻²/year) for each LULC type.

LULC Type	2001	2005	2010	2015	2020
Cropland	264.36	296.30	269.27	296.52	340.77
Forest	290.98	353.76	348.01	367.71	448.65
Grassland	213.69	276.78	277.04	295.75	368.51

Note: LULC refers to land use/land cover.

illustrate the spatiotemporal characteristics of NPP changes across different LULC types. It is evident that NPP in farmlands, forests, and grasslands all exhibited significant upward trends, with forest NPP increasing (290.98 to 448.65 g C·m⁻²/year), showing the most notable growth. This indicates that forest restoration and ecological conservation measures in Shandong Province have made important contributions to enhancing NPP.

Through spatiotemporal dynamic analysis, one can gain a holistic perspective on the transformations occurring in Shandong Province’s NPP and its close relationship with LULC changes can be achieved. These findings provide valuable scientific evidence for ecosystem management in Shandong Province and offer key insights for ecological restoration and sustainable development strategies at the national level.

5.3. The Relationship Between NPP and Human Activities

NPP, a crucial indicator of an ecosystem’s carbon sequestration capacity, is intricately connected to human activities. LULC changes, industrialization, agricultural expansion, pollution emissions, and climate change profoundly affect the global dynamics of NPP [61]. Globally, as urbanization accelerates—particularly in developing countries—large areas of natural vegetation have been replaced by impervious surfaces, significantly reducing NPP levels in these regions [30]. While agricultural intensification can increase NPP in the short term, in the long run, issues such as soil degradation and overexploitation of water resources may suppress sustained agricultural productivity, negatively impacting NPP [62]. On a more positive note, policies such as reforestation and ecological restoration have enhanced NPP in many regions [63], particularly in areas with high forest cover. These ecological projects are vital for sustaining the carbon sink function of ecosystems. Overall, LULC change is the most direct factor influencing NPP, while industrialization and climate change indirectly affect NPP by modifying climatic conditions and ecosystem services.

As a densely populated and economically developed province in eastern China, Shandong Province (2000–2020) has undergone extensive urbanization and industrialization, which has had a profound impact on the region’s ecosystems. The results of this study show that NPP in Shandong Province generally exhibited an upward trend from 2001 to 2020 (Figure 4), but with notable spatial differences across regions. Areas of intense urban development, such as Shandong’s eastern coastal regions and the areas surrounding major cities like Jinan and Qingdao, saw a marked decline in NPP. This decline is primarily due to urban expansion, which has led to widespread coverage by impervious surfaces, including roads, buildings, and infrastructure. These surfaces not only reduce natural vegetation cover but also disrupt the water cycle and soil respiration necessary for plant growth. The PLS−SEM model analysis (Figure 11) indicates a significant negative correlation between NPP and factors such as GDP growth, increased nighttime light intensity, and urban expansion, highlighting the pronounced adverse effects of urbanization on NPP in these regions. In contrast, the central and southern mountainous areas of Shandong Province, benefiting from ecological restoration projects, showed an upward trend in NPP. Projects such as reforestation, wetland restoration, and vegetation rehabilitation have improved the ecological environment in these regions, increasing FVC and significantly enhancing NPP values (Table 4). Notably, the implementation of reforestation projects in the Jiaodong Peninsula and the southern Shandong mountainous areas has led to an increase in forest cover, with these locales exhibiting notably elevated NPP augmentation in comparison to other sectors. This trend is also reflected in LULC changes, as NPP in forests and grasslands has steadily increased during the study period, closely linked to the promotion of ecological protection policies. However, the impact of agricultural activities on NPP in certain areas of Shandong Province cannot be overlooked. Although the expansion of agricultural land has increased NPP to some extent, long-term intensive farming has led to soil degradation and water resource shortages, threatening the productive capacity of ecosystems (Figure 8). This study found that in the highly agriculturalized western Shandong Plain, the increase in NPP was limited, which may be due to over-farming and land degradation. Furthermore, agricultural activities have placed additional pressures on regional water resources, exacerbating vegetation stress, particularly during drought years, when NPP exhibited significant fluctuations. As water resource management policies are implemented, balancing the conflict between agricultural production and ecological conservation will be key to further improving regional NPP. The spatiotemporal variations in Shandong’s NPP are driven not only by direct LULC changes but also by the long-term combined impacts of climate change and human activities. Trend analysis and centroid migration results (Figure 6) show that the NPP centroid has gradually shifted eastward, indicating a significant change in vegetation productivity across both time and space. Meanwhile, the effects of climate factors have shown complex dynamic patterns over different periods, with rising temperatures and decreasing precipitation inhibiting vegetation growth, especially in certain areas (Figure 10). The analysis of Table 2 shows that NPP changes in different regions of Shandong Province exhibit significant spatial variability, with marked improvements in the eastern and central areas, while the western and northern regions have suffered varying degrees of degradation.

In conclusion, human activities have significantly influenced the spatiotemporal variations in NPP across Shandong Province. From urban expansion to agricultural development, LULC changes are the primary driving factors influencing the NPP patterns in this region. Although ecological restoration projects have achieved positive results in some areas, rapid urbanization and agricultural activities remain the main pressures contributing to the decline in NPP. By employing the structural equation model (Figure 11), this study assessed the aggregate impact of multiple determinants on NPP, revealing the intricate interplay between human actions and fluctuations in NPP. These findings provide scientific evidence for future LULC planning and ecological protection policies in Shandong Province and offer valuable insights for sustainable development in other regions.

5.4. Limitations and Future Work

This study, leveraging GEE and multi-source data, examined the NPP dynamics in Shandong Province. Trend analysis and Moran’s I were employed to investigate the spatiotemporal variations of NPP, while PLS−SEM was applied to quantify the influencing factor’s significance. These findings provide important scientific evidence for regional ecosystem management and sustainable development policies. However, there are areas that require further improvement. First, while the CASA model provides an effective framework for NPP estimation, some of the parameter settings in the model rely heavily on existing studies and global averages, which may lead to deviations from actual conditions [64]. This is particularly true across different vegetation types and ecosystems, where variations in light use efficiency could result in regional biases in NPP estimates [65]. Future research could improve the model’s applicability and accuracy by incorporating more field observation data and regional calibration. Second, although this study used high temporal resolution RS data, the spatial resolution of some input data (such as land use type data and meteorological data) remains relatively low. This may limit the ability to accurately capture fine-scale changes in local areas. For example, small-scale LULC changes, microclimate conditions, and local ecosystem characteristics may not be fully reflected in the NPP estimates. Future work should aim to incorporate higher-resolution RS data or other data sources to enhance understanding of the spatiotemporal dynamics of NPP at different scales. Although this study provides a comprehensive analysis of the driving factors of NPP, the multivariable interactions within the model require further investigation. In complex socio-ecological systems, the effects of climate change, LULC changes, and human activities on ecosystems are often non-linear [66], manifesting through a network of intricate interaction mechanisms. Future research could further refine the variable settings in the PLS−SEM model by considering more socioeconomic factors and their interactions with ecosystem productivity. This approach enables a more comprehensive assessment of both the direct and indirect impacts of these factors on NPP.

6. Conclusions

The NPP of vegetation, as a measure of net carbon sequestration, can be utilized to evaluate the potential of vegetation growth in offsetting anthropogenic carbon emissions. It is vital for maintaining ecosystem carbon balance, regulating the carbon cycle, and assessing carbon budgets. This study utilized GEE and multi-source data, along with trend analysis, Moran’s I, and PLS−SEM, to investigate the spatiotemporal dynamics and influential factors of NPP in Shandong Province. The findings reveal that:

(1)

From 2001 to 2020, NPP in Shandong Province exhibited an overall upward trend. This growth was particularly notable in Shandong’s central−southern and eastern regions. The higher values NPP showed concentrated in Shandong’s eastern-southern regions, while lower values were observed in western and northern regions.

(2)

The NPP centroid shifted eastward over the 20-year period, reflecting significant increases in vegetation productivity in the eastern regions. This spatial heterogeneity and centroid shift were closely related to human activities, climate conditions, and topographic features.

(3)

The analysis of driving factors revealed that climatic conditions exerted a significant positive influence on NPP, with temperature and precipitation exerting complex dynamic effects on vegetation productivity across different time periods. Human activity factors, such as urban expansion, nighttime light intensity, and GDP, have been identified to exert a substantial negative influence on NPP, with urbanization in particular accelerating vegetation degradation and reducing regional vegetation productivity. Also, the NPP showed spatial heterogeneity, with elevation and aspect contributing positively to NPP, although their effects varied across different areas.

This study elucidates the regulatory mechanisms underlying vegetation productivity, thereby enhancing our understanding of the spatiotemporal dynamics of NPP in Shandong Province. The findings not only contribute to the theoretical framework of ecological processes but also provide actionable insights for ecological restoration and sustainable development efforts in Shandong. Furthermore, these insights may be extrapolated to other regions globally, offering a template for similar ecological and developmental initiatives.