Spatiotemporal Variation of Net Primary Productivity and Its Response to Climate Change and Human Activities in the Yangtze River Delta, China

Abstract

Exploring the temporal and spatial changes, as well as driving factors, of net primary productivity (NPP) of terrestrial ecosystems is essential for maintaining regional carbon balance. This work focuses on the spatiotemporal variation and future trends of NPP and the response mechanisms of NPP to various driving factors. The Theil–Sen estimator, as well as Mann–Kendall and Hurst exponent methods, were used to analyze the spatiotemporal dynamics and future trends of NPP, and geographical detectors and correlation analysis were used to reveal the response of NPP to various driver changes to environmental factors. The results showed that the NPP was generally on an increasing trend in the Yangtze River Delta region from 2000 to 2019, with the average NPP value of 550.17 g C m⁻² a⁻¹, of which 85.90% was the increasing regions and 14.10% was the decreasing regions, showing a significant spatiotemporal heterogeneity characteristic. The trend of future changes in NPP is dominated by an anti-persistence trend in the study area, i.e., the opposite of the past trend. Notably, annual precipitation is the most significant positive driver of NPP; while NPP was negatively correlated with population, meanwhile, different land use/land cover (LULC) also significantly affected the spatial distribution of NPP. Besides, there was a two-factor enhanced interaction between the various drivers on NPP, with the highest interaction occurring between temperature and elevation. Overall, this study provides data support for future regional NPP predictions and ecosystem evaluations.

1. Introduction

Vegetation is widely distributed in the world and plays an irreplaceable role in the global carbon cycle, water cycle, and material–energy transfer [1,2]. In addition, vegetation also plays a crucial part in regulating the climate, conserving soil and water, and maintaining the stability of ecosystems [3]. Vegetation converts solar energy into biomass through photosynthesis [4,5]. Net primary productivity (NPP) quantifies the total amount of organic matter accumulated by vegetation through photosynthesis and amount remaining after deducting autotrophic respiration per unit of time and space, which reflects vegetation growth and ecosystem health [6], and it is a key factor in evaluating the carbon sink and ecological regulation behavior of terrestrial ecosystems [7,8]. The offset of anthropogenic carbon dioxide (CO₂) emissions by increasing terrestrial NPP is considered to be one of the most environment-friendly and economical ways to mitigate climate change and the greenhouse effect [9,10]. Given the significance of NPP in the global carbon budget and the need for carbon neutral targets on the global scale [11], it is critical to explore the spatiotemporal dynamic characteristics of NPP and its response mechanisms to the driving factors [12].

Changes in NPP are influenced by multiple drivers, including environmental and climatic factors, as well as human activities [13,14,15]. Global climate change has significantly affected ecosystem productivity by altering ecosystem processes, including plant photosynthesis, respiration, phenology, and the rate of soil organic carbon decomposition [16,17,18]. Human activities, such as land use/land cover (LULC) change and urbanization, directly change the structure, processes, and functions of terrestrial ecosystems [6,19], which consequently affect terrestrial NPP [20,21,22]. Recently, many studies have been conducted on understanding the response of NPP to different factors at different temporal and spatial scales, ranging from inter-monthly or inter-annual to millennial, as well as from landscape and regional to global [6,14,22,23,24,25], with methods such as correlation analysis [26], regression [27], principal component analysis [28], sensitivity analysis [29], etc. However, there are still uncertainties in the relative impacts of climate change and human activities on NPP, although several studies have focused on these issues and attempted to explain the response of NPP to climate change and anthropogenic perturbations [7,13]. Previous studies have taken two main approaches to assess the relative contributions of climate change and human activities to terrestrial ecosystems NPP [30,31,32]. The first approach uses regression models (pixel-specific linear regression [33], partial least squares regression (PLSR) model [34], multiple linear regression [35], etc.) to assess the relative contribution of each factor [36]. Another approach is to compare the residuals between observed NPP and predicted NPP by modeling methods, such as land surface modeling [37]. However, quantitative assessment of the relative impacts of climate change and human activities on spatial-scale changes in NPP remains a challenge for NPP research in terrestrial ecosystems due to the lack of reasonable assessment methodologies [38]. In the face of changing climate and on their implications in forest productivity, temperature and drought have great effects on forest growth under different altitudes and phenological conditions, elevation gradient is an important indicator of ecological process [39,40]. Overall, the interpretation of the combined effects of multiple factors on NPP, especially the interaction between these factors, has yet remained unclear [12].

The Yangtze River Delta is located in the central-eastern part of China. It is one of the most dynamic regions in the global economy and one of the most urbanized regions in China [41]. Over the past decades, the Yangtze River Delta region has experienced rapid industrialization and urbanization, while rapid economic development and intensive land development have resulted in a dramatic decline in regional ecosystem quality and biodiversity [42]. Moreover, the Yangtze River Delta is considered to be sensitive to climate change due to the influence of the Asian monsoon [43]. For example, Wu et al. [12] used spatiotemporal analysis to determine the relationship between NPP, nighttime light urbanization index values, and climatic factors from pixel to regional scales. The results revealed that the urbanization factor is the main driving force for NPP change in high-speed urbanization areas, and the factor accounted for 47% of the variations. However, in the forest and farm regions, the NPP variation was mainly controlled by climate change and other factors. Therefore, exploring the spatiotemporal dynamic characteristics of long-term NPP, in this region, can provide insight into the effects of human activities and climate change on ecosystem structure and function and, thus, provide decision-making data support for regional urban planning and ecological construction.

There are many papers which have investigated the change in the carbon cycle in the research region. However, up to now, there are few published research works which focus on the spatiotemporal variation and future change trends of NPP and the response mechanisms of NPP to various driving factors concerned with the research region. Therefore, we make two hypotheses in the manuscript: (1) The future change trends of NPP in the research region are the same as the past; (2) among the several driving factors impacting NPP, human activities are more important than climate factors.

2. Materials and Methods

2.1. Study Area

As an alluvial plain formed before the Yangtze River enters the sea, the Yangtze River Delta region is located in the lower reaches of the Yangtze River in China, bordering the Yellow Sea and the East China Sea, and at the confluence of the rivers and the sea. The region includes a total of 41 cities in Shanghai, Jiangsu, Zhejiang, and Anhui provinces, with a total area of 358,000 km² (Figure 1). “Construction land” refers to the land for building buildings and structures [44]. By 2020, the Yangtze River Delta region had a population of 227 million and a gross domestic product (GDP) of 24.5 trillion RMB, and the urbanization rate of the resident population exceeded 60%. With less than 4% of the country’s land area, the region creates nearly 1/4 of China’s total economic output and 1/3 of its total imports and exports [45].

The topography of the Yangtze River Delta is low and flat, with elevations, basically, below 10 m above sea level. Most of the lands are lower plains, with some hills and isolated mountains scattered in the southern part of the study area. The research region is located in the distribution area of a subtropical monsoon climate, which is influenced by the Asian monsoon and has a humid and mild climate with four distinct seasons and sufficient rainfall, with an average annual temperature of 16.90 °C, and with an annual precipitation of about 1235 mm (Figure 2). The regional vegetation cover types mainly include subtropical evergreen broad-leaved forests, mixed deciduous broad-leaved forests, coniferous forests, etc.

2.2. Data Source and Preprocessing

The data sources of this study mainly include NPP data, human activities data, and meteorological data. The version of the NPP data is MOD17A3HGF v006 and was obtained from the Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov, accessed on 20 December 2021), with the spatial resolution of 500 m [46]. In addition, this data product contains a data quality control file (NPP_QC), which indicates the quality reliability of NPP products in different regions. We conducted statistical analysis based on the NPP_QC data from 2000–2019, and the results showed that the multi-year average reliability of medium and high-level NPP data quality reached 90.32% in the Yangtze River Delta region from 2000 to 2019 [47]. The LULC data were obtained from the Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center (LAADS DAAC, https://ladsweb.modaps.eosdis.nasa.gov/, accessed on 2 March 2020), for MCD12Q1 v006 data, from 2001 to 2019. The population density [48] and GDP data (all with 1 km spatial resolution) were provided by the Data Center for Resources and Environmental Sciences of the Chinese Academy of Sciences (https://www.resdc.cn, accessed on 5 March 2022). We downloaded the digital elevation model (DEM) data [49] (GDEMV3 30M) with 30 m spatial resolution from the Geospatial Data Cloud (http://www.gscloud.cn, accessed on 5 March 2022), and they were processed using ArcGIS10.7 software to obtain slope and aspect data. The data sources for the spatial data used in this study are shown in

Table 1. The data sources used in this study.

Category	Data	Year	Original Resolution
Net primary productivity data	MOD17A3HGF v006	2000–2019	500 m
Land use/land cover data	MCD12Q1 v006	2001–2019	500 m
Population data	Population density	2000, 2005, 2010, 2015, 2019	1000 m
Gross domestic product data	GDP
Digital elevation model data	GDEM v003	2019	30 m
Meteorological data	Temperature, Precipitation, Sunshine hours data	2000–2019	500 m

.

The meteorological data are obtained from the monthly dataset of China’s terrestrial climate information provided by the National Meteorological Information Center (http://data.cma.cn, accessed on 10 February 2022), which mainly includes monthly average temperature, total precipitation, and total sunshine hours from 94 meteorological stations in and around the Yangtze River Delta region. The data period was 2000–2019, and the annual average data of climate factors at each station were obtained by aggregating monthly data, and the thin plate spline method was used to spatially interpolate the meteorological data within the region. All the above data were converted to raster data with the same projection coordinate system and 500m spatial resolution in ArcGIS 10.7 after a series of data pre-processing, including mosaic, projection transformation, resampling, and clipping. In addition, in order to meet the needs of the geographical detector input variable types (i.e., categorical variables), each driving factor was discretized and classified into nine categories using the natural breaks method [50], as shown in

Table 2. The spatial driving factors of NPP in this study.

Category	Environmental Driving Factors			Climatic Driving Factors			Human Activities
	DEM/m	Slope/°	Aspect	TEM/°C	PRE/mm	TSH/h	GDP	POP	LULC
1	−56~66	0~0.60	Flat	8.90~11.04	766~959	1450~1598	92~63,114	45~316	Forest land
2	66~182	0.60~1.80	North	11.04~12.81	959~1102	1598~1669	63,114~126,137	316~984	Grassland
3	182~311	1.80~3.26	Northeast	12.81~14.21	1102~1232	1669~1731	126,137~189,159	984~2320	Construction Land
4	311~441	3.26~4.71	East	14.21~15.24	1232~1368	1731~1799	189,159~252,182	2320~4591	Cultivated land
5	441~577	4.71~6.25	Southeast	15.24~15.94	1368~1499	1799~1873	252,182~315,204	4591~8198	Water
6	577~726	6.25~7.97	South	15.94~16.51	1499~1610	1873~1950	315,204~378,227	8198~13,274	Other
7	726~894	7.97~9.94	Southwest	16.51~17.21	1610~1728	1950~2024	378,227~441,249	13,274~18,618	/
8	894~1101	9.94~12.51	West	17.21~18.12	1728~1939	2024~2098	441,249~504,272	18,618~23,962	/
9	1102~1593	12.51~21.84	Northwest	18.12~19.39	1939~2343	2098~2237	504,272~567,294	23,962~34,115	/

Notably, “TEM” represents the annual average temperature, “PRE” represents the total annual precipitation, “TSH” represents the total annual sunshine hours, and “POP” represents the number of populations.

. The natural breaks, also called the Jenks natural breaks classification method, is a data clustering method designed to determine the best arrangement of values into different classes. This is done by seeking to minimize each class’s average deviation from the mean value of the class, while maximizing each class’s deviation from the means of the other classes. In other words, the method seeks to reduce the variance within classes and maximize the variance between classes. Therefore, it is very applicable to the categorization of the various driving factors.

2.3. Methods

2.3.1. Trend Analysis

The Theil–Sen estimator method is a calculation method applicable to trend analysis for long-term series data, which is computationally efficient, insensitive to measurement errors and niche data, and can be combined with the Mann–Kendall trend analysis method to further perform trend significance tests for long-term series data [51]. Its calculation formula is as follows:

$$\beta =median \left(\frac{NP{P}_j-NP{P}_i}{j-i}\right)$$

(1)

where: median () represents taking the median value; NPP_i and NPP_j are the corresponding data (NPP) values in years i and j, respectively, and j > i; when β = 0, NPP has no change, when β > 0, NPP increases, and vice versa.

The Mann–Kendall trend test method is widely used for trend analysis of NPP time series, as it does not require the measurements to obey a normal distribution and is not affected by missing values and outliers when performing trend significance tests [52]. The process is as follows: for the sequence X_t = x₁, x₂,..., x_n, first determine the relationship between the magnitude of x_i and x_j (set to S) for all pairs of values (x_i, x_j, j > i). Hypothesis the data in the H₀ series are randomly aligned, i.e., there is no significant trend; there is an upward or downward trend in the H₁ series. The test statistic S is calculated as:

$$S={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^{n}sgn\left(x_j-x_i\right)$$

(2)

$$sgn\left(x_i-x_j\right)=\left\{\begin{array}{c}1, x_j-x_i>0\\ 0, x_j-x_i=0\\ -1, x_j-x_i<0\end{array}\right.$$

(3)

where: x_i and x_j are the corresponding data (NPP) values of the i-th and j-th years, respectively, and j > i; n is the length of the data set. Then, there is

$$Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{Var\left(S\right)}}, S>0\\ 0, S=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}}, S<0\end{array}\right.$$

(4)

where: Z is the test statistic; Var(S) is the variance. The critical value Z₁−α/2 was found in the normal distribution table at a given significance level. When |Z| ≤ Z₁−α/2, the original hypothesis is accepted, i.e., the trend is not significant and vice versa. In addition, we used the piecewise regression model to identify the turning point of interannual variation in NPP [6].

2.3.2. Future Change Trend Analysis Based on Pixels

The Hurst exponent [53], based on the rescaled range (R/S) analysis method, is an important indicator for quantitatively describing the long-term dependence of time series data, and its basic principle [54,55,56] is that, with a time series ξ(t) (t = 1, 2…), for any positive integer $\tau$≥ 1, there is a mean series:

$$<\xi {>}_t=\frac{1}{\tau }{\displaystyle \sum }_{t=1}^{\tau }\xi \left(t\right) \left(\tau =1,2,\cdots \right)$$

(5)

The cumulative deviation X (t) is thus obtained as:

$$X\left(t,\tau \right)={\displaystyle \sum }_{t=1}^t\left[\xi \left(t\right)-<\xi {>}_t\right] \left(1\le t\le \tau \right)$$

(6)

The polar difference R is defined as:

$$R\left(\tau \right)=\underset{1\le t\le \tau }{\max }X\left(t,\tau \right)-\underset{1\le t\le \tau }{\min }X\left(t,\tau \right) \left(\tau =1,2,\cdots \right)$$

(7)

The standard deviation S is defined as:

$$S\left(\tau \right)={\left[\frac{1}{\tau }{\displaystyle \sum }_{t=1}^{\tau }{\left(\xi \left(t\right)-<\xi {>}_{\tau }\right)}^2\right]}^{\frac{1}{2}} \left(\tau =1,2,\cdots \right)$$

(8)

R, S, and τ satisfy the general relationship:

$$R\left(\tau \right)/S\left(\tau \right)=c·{\tau }^H$$

(9)

Least squares fitting:

$$\log {\left(R/S\right)}_{\tau }=\log c+H·\log \tau$$

(10)

where: c is a constant; R(τ)/S(τ) is the rescaled range; H is the Hurst exponent. The Hurst exponent indicates the range of the cumulative deviation from the mean of the time series over time. The Hurst exponent of the time series is between 0 and 1. There are three main scenarios: (1) 0 < H < 0.5, indicating that the future trend of the vegetation NPP is anti-persistence, i.e., opposite of the past trend; (2) 0.5 < H < 1, indicating that the future trend of the vegetation NPP is persistence, i.e., consistent with the past trend; (3) H = 0.5, and the time sequence is random.

2.3.3. Correlation Analysis

In this study, the Pearson correlation coefficient, based on pixels, was used to assess the response relationship between the annual NPP of vegetation to different driving factors, which was calculated as follows:

$$r_{x,y}=\frac{{\sum }_{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(11)

where, r_x,y is the x, y correlation coefficient; x_i and y_i represent the values of variables x and y in period i, respectively; $\overline{x}$ and $\overline{y}$represent the average values of variables x and y; the value of the correlation coefficient ranges from negative to positive, where a positive value represents positive correlation between variables, a negative value represents negative correlation between variables, and the larger the absolute value of correlation coefficient, the stronger the correlation between variables.

2.3.4. Geographical Detectors

Geographical detectors are statistical methods to detect spatial heterogeneity and reveal its underlying driving forces, which consist of four modules: factor detection, risk detection, interaction detection, and ecological detection [57]. Geographical detectors enable the assessment of the extent of impact of multiple driving factors on response variables and the magnitude of interactions between factors. Among them, the factor detection and interaction detection modules are widely used in the analysis of driving factors of geographic phenomena and processes, and the factor detector detects the spatial heterogeneity of the dependent variable and the extent to which the independent variable explains the spatial heterogeneity of the dependent variable, as measured by the q-value. The expression is:

$$q=1-\frac{{\sum }_{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }^2}$$

(12)

where: h = 1,..., L is the stratification of variable Y or factor X, i.e., classification or partitioning; N_h and N are the number of cells in stratum h and full area, respectively; ${\sigma }_h^2$ and ${\sigma }^2$ are the variances of Y values in stratum h and full area, respectively. q-values are between 0 and 1, and larger values indicate that the explanatory factor has a stronger effect on the study variable the stronger the driving force.

The greatest advantage of the geographical detector, in relation to other statistical methods, is its possession of interaction detectors that can be used to identify interactions between different risk factors, i.e., to assess whether two driving factors, together, increase or decrease the explanatory power of the dependent variable or whether the effects of these factors on the dependent variable are independent of each other [58].

The risk detector is used to detect the suitable range or type of impact of different environmental factors on NPP, which are tested by t statistics:

$$t=\frac{{\overline{Y}}_{h=1}-{\overline{Y}}_{h=2}}{{\left[\frac{Var\left(Y_{h=1}\right)}{n_{h=1}}+\frac{Var\left(Y_{h=2}\right)}{n_{h=2}}\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}$$

(13)

where: ${\overline{Y}}_h$ is the mean value of NPP in subregion h; $n_h$ is the number of samples in subregion h; $Var$ is variance.

3. Results

3.1. Spatiotemporal Variation Characteristics of NPP in Yangtze River Delta

3.1.1. Spatial Distribution Pattern of NPP

The average NPP was 550.17 g C m⁻² a⁻¹ in the Yangtze River Delta from 2000 to 2019, with NPP values ranging from 68.21 to 1433.10 g C m⁻² a⁻¹ (Figure 3). Among them, the percentage of regions with average NPP values greater than 400 g C m⁻² a⁻¹ was 93.04%, and the percentage of more than 600 g C m⁻² a⁻¹ was 30.93%. Generally, the NPP in the research region shows the spatial distribution characteristics as low in the northwest and high in the southeast, and the high value areas were mainly distributed in the mountainous and hilly areas with relatively high elevation (Figure 3), such as southern Anhui Province and south-central Zhejiang Province. In particular, the NPP of Zhejiang Province was significantly higher than that of other provinces and cities, which may be due to the high forest cover and abundant vegetation species resources in Zhejiang Province. In addition, there is less impact on vegetation from human activities in the region, and vegetation enables natural growth. Low value NPP areas were mainly located in north-central Anhui Province, southern and northern Jiangsu Province, and northern Zhejiang Province, including Xuzhou city, Huaibei city, Luan city, Hefei city, Suzhou city, and the border area between Hangzhou city, Jiaxing city, and Huzhou city. These areas are flat and low in elevation, and the main land types are cultivated and construction land (Figure 1), which explains the low NPP values in this area.

3.1.2. Temporal Variation Characteristics of NPP

As shown in Figure 4, the average NPP showed an oscillating upward trend, at a rate of 3.3 g C m⁻² a⁻¹, in the research region during 2000–2019, with a peak value of 609.84 g C m⁻² a⁻¹ in 2014, followed by a decrease and then an increase. According to the piecewise regression model [6], the year of 2009 was extracted as the turning point of the NPP trend. From 2000 to 2009, the average NPP increased by 46.84 g C m⁻² a⁻¹ in the Yangtze River Delta, and the annual growth rate of NPP was 4.68 g C m⁻² a⁻¹. From 2010–2019, the growth rate of NPP in the Yangtze River Delta slowed down to 3.55 g C m⁻² a⁻¹, with an average NPP growth of 35.5 g C m⁻² a⁻¹ over the 10-year period.

In addition, the results of the Theil–Sen estimator for NPP, during the period of 2000–2019, indicate that the area of positive β values was significantly larger than the area of negative values. Specifically, 85.90% of the regions had increasing trends and 14.10% of the regions had decreasing trends, and the mean value of β was 1.63, which indicated that the NPP in the research region had been in an overall increasing trend in the past 20 years. We further performed the Mann–Kendall test on the NPP trends, and the results showed that there were significant spatial heterogeneity characteristics of NPP trends in the Yangtze River Delta (Figure 5). Among them, the percent of significant increase in the region accounted for 47.25%, and the percent of significant decrease was only 5.16%; no significant increase region accounted for 35.08%, and no significant decrease accounted for 12.51%. The percentage of regions with no change was extremely small and negligible. We observed that the significance of NPP in most of the research regions was in increasing trends (including significant increase and non-significant increase), and the significant decrease was sporadically distributed in the southern part of Zhejiang Province (e.g., Lishui city, Wenzhou city, Taizhou city), southeastern Jiangsu Province (e.g., Suzhou city, Nan-tong city), and Shanghai city.

3.1.3. Future Change Trends of NPP

We analyzed the future change trend of NPP in the research region using the Hurst exponent, at the pixel level, based on the NPP data from 2000–2019. The Hurst exponent was further divided into five classes according to the value range (Figure 6).

The future trend of NPP, in the research region, was 30.83% in the persistent region and 69.17% in the anti-persistent region, in which the proportions of weak anti-persistence and weak persistence were 68.89% and 30.41%, respectively, which indicated that the future trend of the research region was mainly in the anti-persistence state, i.e., the opposite of the past trend. From the spatial distribution of NPP future trends, weak anti-persistence was widely distributed in most regions, while weak persistence was concentrated in northern Anhui and urban clusters near the Yangtze River valley.

3.2. Relationship between Various Driving Factors and NPP

The driving factors of spatiotemporal variation of NPP in the research region are analyzed using the geographic detectors, and the indicator q value of the factor detectors can be used to measure the intensity of each driver’s influence on NPP (Figure 7). As can be seen from Figure 6, annual precipitation is the most important driving factor affecting NPP in the research region, followed by elevation, slope, annual average temperature, and LULC. To explore the response of NPP to different driving factors, we perform further detailed analysis from the perspective of different types of drivers.

3.2.1. Relationship between Environmental Driving Factors and NPP Change

Risk detectors can reflect the appropriate type or range of NPP distribution of various factors, providing scientific basis for ecological protection and restoration. As the effect of the Aspect factor on NPP was low (q-value of 0.002) (Figure 7), it was not considered in the analysis of the impact of environmental driving factors on NPP changes (same for GDP). There was a similar change pattern between elevation and slope classification and NPP data, i.e., the value of NPP increases with the increase in elevation or slope (Figure 8), which could be due to the fact that the slope factor was extracted from the elevation factor. Specifically, after the classification of elevation or slope reached class 4, the increasing trend of NPP slowed down significantly until class 8, and at class 9, the NPP corresponding to elevation maintained an increase, while the NPP corresponding to slope decreased.

3.2.2. Relationship between Climate Driving Factors and NPP Change

According to the results of the Pearson correlation analysis of NPP, NPP was positively correlated with annual average temperature (p < 0.001) and annual precipitation (p < 0.001), and it was negatively correlated with sunshine hours (p < 0.001) in the north-central Yangtze River Delta (Figure 9). The NPP was negatively correlated with annual precipitation and positively correlated with sunshine hours in the southern Yangtze River Delta, which was partially positively and partially negatively correlated with annual average temperature.

In particular, the average Pearson correlation coefficient between NPP and annual average temperature in the research region was −0.0037, the area of the positively correlated area accounts for 51.12%, the area of the negatively correlated area accounts for 48.88%, and the ratio of positively and negatively correlated area was close to 1:1 (Figure 9a). The average Pearson correlation coefficient between NPP and annual precipitation was 0.18, and the area of positively correlated areas accounted for 78.88% and 21.12% were negatively correlated areas, where annual precipitation generally contributed to the growth of NPP (Figure 9b). The average Pearson correlation coefficient between NPP and sunshine hours was −0.06, with 37.27% of the area positively correlated and 62.73% of the area negatively correlated (Figure 9c). Overall, NPP in the research region showed a strong positive correlation with annual precipitation, and annual precipitation contributed to an increase in NPP, except for the southern region of the study area.

3.2.3. Relationship between Human Activities and NPP Change

The average Pearson correlation coefficient between NPP and population in the Yangtze River Delta was −0.0096, the positively correlated area accounts for 35.22%, (e.g., Xuancheng city, Huzhou city, Hangzhou city, and Jinhua city, etc.), and the negatively correlated area accounts for 64.78%, mainly in the northwestern and southern parts of the research region (Figure 10).

Different land cover types correspond to different land ecosystem structures, their NPP are also different, and the variability of land cover types and their spatial distribution directly determines the spatial distribution characteristics of NPP. As illustrated in Figure 11, considering the NPP corresponding to each land type during 2001–2019, forest land NPP (655.51 g C m⁻² a⁻¹) was the highest, followed by grassland NPP (600.12 g C m⁻² a⁻¹) and cropland NPP (487.51 g C m⁻² a⁻¹). Despite the large differences in NPP among different land use types, the NPP trends of forest land, grassland, and cultivated land were similar, with all showing an oscillating upward trend in which cultivated land had the fastest rate of increase, followed by forest land and grassland.

3.2.4. Interaction of Various Driving Factors on the Effect of NPP

The changes in NPP often result from a combination of various factors, and we further explored the interaction of the effects of different drivers on NPP in the research region (Figure 12). The interaction between two factors was greater than that of a single factor. In particular, the interaction of annual average temperature and elevation had the highest q-value of 0.607, which was greater than the single factor q-value of annual precipitation, which indicated that the interaction of annual average temperature and elevation exceeded the effect of annual precipitation on NPP. In addition, the q-values between annual precipitation and other factors were all greater than 0.516 (q-value of annual precipitation as a single factor), which implied that the synergistic effect of other factors with annual precipitation increased the explanatory power of the effect of annual precipitation on NPP. This further validates that annual precipitation is the most significant driving factor for the Yangtze River Delta.

4. Discussion

In this study, we explored the spatiotemporal variation of NPP in the Yangtze River Delta, from 2000 to 2019, and the relationship between various driving factors and NPP. The results showed that NPP increased with elevation and slope in the research region. In general, the forest land use type has a higher NPP [52], and it is mainly found in higher elevations in mountainous and hilly areas (Figure 1 and Figure 3). Under a certain elevation (e.g., below 1593 m elevation in this study area), NPP increased rapidly with elevation and then maintained a steady increase (Figure 8). There is more interest in the impact of climate variables and human activities on NPP than other environmental variables. Therefore, this paper focuses on the effects of climate change and human activities on NPP and the future change trend of NPP in the region.

4.1. Impact of Climatic Factors on NPP

Due to its special geographical location, the research region is vulnerable and sensitive to climate change, as it is affected by both sea–land thermal differences and seasonal variations in circulation [12,59,60]. Some studies have shown that the growing season in the research region has experienced a dramatic warming and drying trend in the past 20 years [60]. However, the area of the research region is very large (358,000 km², Figure 1), and the difference in elevation and climate are relatively significant within the region, e.g., the spatial distribution of temperature, annual precipitation, and sunshine hours are significantly different. This is reflected in the significant spatial heterogeneity of Pearson correlation coefficients of NPP with mean annual temperature, annual precipitation, and sunshine hours.

Previous studies have shown that NPP increases with temperature in many regions [14,18]. However, this finding does not completely hold true in the research region. The Pearson correlation analysis results for annual average temperature and NPP showed that more than 48.88% of the regions exhibited a negative correlation between NPP and annual average temperature, and about 51.12% of the regions showed a positive correlation between NPP and annual average temperature (Figure 9a). The results of these analyses show that the warming of nearly half of the research region does not directly contribute to the increase in NPP. This was likely due to the fact that the research is in the southern region of China, which belongs to the subtropical monsoon climate, so the main limiting factor of vegetation growth is not temperature [12]. However, increased temperature may increase soil evapotranspiration, cause soil moisture stress, reduce soil water content, and thus affect the vegetation growth process [61,62,63].

Annual precipitation was the most dominant positive factor for NPP changes (Figure 9b), with more than 78% of the area positively correlated with annual precipitation (Figure 8b), which indicated that an increase in annual precipitation, in most areas of the study region, would contribute to an increase in NPP. Our results are consistent with previous studies [12,60]. In addition, the areas negatively correlated with annual precipitation (about 21.12%) are mainly located in the higher elevation mountains and hills in the southern part of the research, which may be due to the higher vegetation cover in this region, and the vegetation has the functions of water conservation, precipitation retention, runoff regulation, and evapotranspiration reduction [64]. Thus, the increase in precipitation does not result in the increase in NPP in this region.

Most plants have their own biorhythmic clocks, which can sense changes in the external light environment. Suitable sunshine hours are the basis for effective photoenergy absorption, transfer, and conversion, which affect the growth of plants [65]. Photosynthesis of vegetation increased with the increase in sunshine hours. In contrast to annual precipitation, vegetation in the southern part of the study area was positively correlated with sunshine hours (Figure 9c), and increasing sunshine hours in this region could increase NPP. Moreover, NPP was negatively correlated with sunshine hours in the northern part of the study area, and this phenomenon was also observed in the Pearl River Delta region [6]. The main vegetated land type in the northern part of the study area is cultivated land (Figure 1), where the main crops are wheat and corn, and the sunshine hours in the area may have met their growth requirements [66].

4.2. Impact of Human Activities on NPP

In the past 20 years, with the increase in the range and intensity of human activities, anthropogenic disturbance has become a non-negligible driving force affecting NPP in the research region [6,12]. In the research region, the most significant forms of natural transformation by human activities are changes in land cover. From 2000 to 2019, the NPP of different vegetation land use types in the research region showed an oscillating upward trend (Figure 11). Rapid urban expansion and economic development have caused tremendous changes in the LULC and the ecological environment surrounding the city [6]. From 2000 to 2019, the NPP of different vegetation types showed an oscillating upward trend in the Yangtze River Delta, which was attributed to the implementation of the Grain for Green Project, Yangtze River Protected Forest Project, and Ecological Restoration Project [13,67]. Therefore, increasing the area of urban green space and rational planning of urban structure could contribute to the increase in NPP. Nevertheless, some studies have shown that the Yangtze River Delta region has lost 0.14 Tg C of NPP per year due to urbanization [12].

In addition, the impact of urbanization on NPP not only has an impact on land use type changes but also on the rapid growth of urban population. Over the past several decades, China’s urbanization rate has increased dramatically from 17.92% to 59.58% [68], and the rapidly increasing urbanization rate may represent a huge population growth. Previous studies have confirmed that the gross primary productivity (GPP) and enhanced vegetation index (EVI) both showed a significant decreasing trend with urban expansion [69]. The dramatic LULC changes caused by rapid urbanization have led to a decline in NPP, which is mainly due to the replacement of cropland and forests by urban land (i.e., construction land) [70]. The NPP was negatively correlated with population in more than 64% of the regions (Figure 9), which indicated that population growth in the region is one of the major factors resulting in the decrease in NPP. Previous studies have shown there is a negative correlation between NPP, economic expansion, and urbanization. As urban population increases, cities have been expanded to surrounding green spaces (including forests, grasslands, croplands, etc.) to accommodate more people, which, in turn, results in the decrease in NPP [71]. Wu et al. calculated the NPP in Hubei province between 2001 and 2012 and studied the spatial and temporal dynamics of the NPP and its relationship with urbanization indicators and found that there was a negative correlation between NPP, population, and GDP. In general, the highest values of negative correlation coefficients were found in urban areas, and the lowest values were found in woodlands. The negative correlation coefficients between population and annual NPP became larger as population density increased [72]. This might explain the negative correlation between population density and NPP in the research region. In recent years, the research region has proposed the “Yangtze River Delta Regional Integrated Development Plan”, which may help improve cross-basin and cross-regional ecological environmental quality and further enhance regional NPP as the policy is implemented [73]. Therefore, controlling population growth and planning urban development appropriately will contribute to improving the NPP of the study area.

4.3. Future Change Trends of NPP in Yangtze River Delta

In the past 20 years, the NPP was generally on an upward trend in the Yangtze River Delta, and most of the regional NPP values were higher than the national average NPP [6]. However, the uncertain implications of future climate change [2] and human activities [13] could result in changes in the trend of NPP in the Yangtze River Delta. Our results show that 68.89% and 30.41% of the NPP future trends in the research region were weak in anti-persistence and weak in persistence, respectively (Figure 6), indicating that the future trends in the research region were mainly in the anti-persistence state, i.e., opposite to the past trends.

To further explore the persistence of future NPP trends in the research region, we spatially overlaid the Hurst exponent and Mann–Kendall trend analysis of NPP layers to obtain eight future NPP change scenarios (Figure 13). As Figure 13 illustrates, the highest percentage of the future trend of the NPP, in the research region, is a significant increase in weak anti-persistence with 32.74%, followed by no significant increase-weak anti-persistence (26.00%) and significant increase-in weak persistence (14.51%). The no significant increase and the significant increase trends were dominated by weak anti-persistence in most areas, which indicated that the future trend of NPP might be dominated by a decrease in the study area. This opposite pattern of NPP change from the past may be influenced by a combination of climate change and human activities [13,14]. The research region has sufficient water and heat conditions, as ~20–22 degrees Celsius is the optimum growth temperature for vegetation in the region [74], and increased annual precipitation contributes to the increase in NPP in most areas. However, the impact of anthropogenic disturbances on NPP should not be underestimated as well. In the urbanized area of East China, anthropogenic disturbances caused a decrease in the average NPP from 818 C m⁻² a⁻¹ in 1991 to 699 C m⁻² a⁻¹ in 2002 [75]. Urbanization generally decreases NPP; however, NPP can be increased by introducing high productivity plants, increasing urban green areas, and reducing temperature, water, and nutrient limitations [76].

4.4. Strengths and Uncertainties

This work assesses the spatiotemporal variation characteristics and future trends of NPP in the research region and the impact of various drivers on NPP in 20 years based on geographical detectors. Unfortunately, the geographical detector can only test the coupling of the spatial distribution of two variables and is not able to detect the coupling between three or more driving factors [57,58]. Furthermore, only nine driving factors were considered in this paper: Indeed, numerous climatic factors and human activities also affect the vegetation growth process [77,78], and more determinants should be included in future studies, e.g., solar radiation, soil type, and ecological engineering projects [79,80].

The time lag effect of NPP on various driving factors is also not considered in this paper, which should be further studied on different time and spatial scales in the future [14,77,78]. Besides, only one NPP data product was used in this paper, and further work should include more types of NPP to reduce the uncertainty in assessing the spatiotemporal variability of NPP and driving force analysis [81,82]. Despite the above limitations, our study reveals the spatiotemporal variation characteristics of NPP, future change trends, and the response of NPP to various driving factors in the research region during 2000–2019. The results of these analyses provide a reference for coping with global climate change and regional socio-economic development.

5. Conclusions

The NPP was generally on an increasing trend in the research region from 2000 to 2019, which showed significant spatiotemporal heterogeneity characteristics, i.e., low in the northwest and high in the southeast. The future NPP change trend is dominated by anti-persistence, i.e., the opposite of the past trend.

Among the several driving factors impacting NPP, Annual precipitation was the most significant positive climate driver of NPP. Land use/cover had a significant impact on NPP changes, and urban expansion and population increase resulted in the decrease in NPP. In addition, there was a two-factor enhanced interaction between the impact of various drivers on NPP, with the highest interaction between annual average temperature and elevation. In summary, the NPP was impacted by a combination of climate change and human activities, and it could be increased in the study area by controlling population growth, appropriately planning urban development, and increasing the area of green vegetation.

There are still many other factors to be considered in future studies, such as solar radiation, soil type, vapor pressure deficit (VPD), and ecological engineering. In addition, the geographical detectors can only test the coupling of the spatial distribution of two variables (i.e., the relationship between two variables), and they are unable to detect the coupling between three or more driving factors. Therefore, other analytical approaches (e.g., structural equation modeling) are needed to further explore and assess the response mechanisms of NPP to various driving factors so that the understanding of the links between NPP and climate change can be strengthened in the future.